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ABSTRACT 



Finite Temperature QCD with Domain Wall Fermions 
GEORGE TAMMINGA FLEMING 

Domain wall fermions are a new lattice fermion formulation which preserves the 
full chiral symmetry of the continuum at finite lattice spacing, up to terms exponen- 
tially small in an extra parameter. We discuss the main features of the formulation 
and its application to study of QCD with two light fermions of equal mass. We also 
present numerical studies of the two flavor QCD thermodynamics with aT — 1/4. 
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Introductory remarks 



For the last fifty years, quantum field theories have been amazingly successful at 
predicting observable physical phenomena involving elementary particles. A class of 
quantum field theories called gauge theories, because each contains a special inter- 
nal symmetry called gauge symmetry, accurately describes the interactions, or forces, 
among elementary particles. Three of the four known fundamental forces {strong, 
weak, and electromagnetic) are best described by gauge theories. Furthermore, nearly 
all promising new theories are based on the principle of invariance under local sym- 
metry transformations. 

Quantum Chromodynamics, or QCD, is the gauge theory of the strong force re- 
sponsible for binding protons and neutrons together to form nuclei. As a practical 
matter, although QCD is the underlying microscopic theory responsible for nuclear 
physics, theorists do not use QCD to study the formation of nuclei. The calculation is 
possible in principle, but the large number of degrees of freedom render it intractable. 
However, when particles interact at high energies, the strength of the strong interac- 
tion becomes paradoxically weak through the phenomenon called asymptotic freedom. 
Only in this perturbative regime can expansions of QCD in Feynman diagrams be 
used to successfully describe observed physical processes. 

Another approach to understanding QCD employs the formal mathematical equiv- 
alence between statistical mechanics and quantum field theories expressed in the 
Euclidean path integral formalism. Systems described by statistical mechanics also 
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contain vast numbers of degrees of freedom yet numerical and analytic techniques 
have been developed to handle them. When QCD is expressed in the Euchdean path 
integral formalism, the analogy with statistical mechanics becomes manifest with the 
QCD Lagrangian as the Boltzmann weight. Unfortunately, the path integral over a 
continuous spacetime is divergent unless the integrals are properly regularized. 

Lattice QCD provides this regularization scheme by discretizing the Euclidean 
spacetime manifold while exactly preserving the internal gauge symmetry. Unfortu- 
nately, the continuous Lorentz symmetries of translation and rotation are explicitly 
broken down to discrete subgroups in the process. The goal is to study QCD on as 
fine a lattice as possible to understand the continuum limit of infinitesimal lattice 
spacing where Lorentz symmetries are restored. Due to asymptotic freedom, the ra- 
tios of physical observables will become independent of the lattice spacing when the 
spacing is small enough. For the lattice QCD theorist, accomphshing a few steps to- 
ward the continuum limit is an immense undertaking involving large supercomputers. 
The reward for this effort is the ability to study non-perturbative regimes of QCD. 

QCD at finite temperature is a non-perturbative regime where lattice QCD has 
emerged as an ideal tool for direct study of the thermodynamics. Theorists have 
long speculated that QCD at sufficiently high temperature or density undergoes a 
deconfining phase transition where quarks and gluons are no longer confined within 
meson and baryon bound states. In the limits of either infinite or zero mass quarks, 
the phase transition is characterized by the behavior of a physical observable called 
an order parameter. In nature, though all the quarks have finite mass, the up and 
down quarks are much lighter, as evidenced by the light pion masses. So, it is likely 
that QCD in our world resembles QCD with two massless quarks. When the quarks 
are massless, the deconfinement order parameter is associated with the spontaneous 
breakdown of chiral symmetry, the mechanism that generates massless pions. 

In this work, we use a recently developed fermion formulation, domain wall 
fermions, to study aspects of QCD that are strongly infiuenced by chiral symmetry. 
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Using fermion formulations that were available in the past, the full chiral symmetry 
of QCD is broken at finite lattice spacing along with the Lorentz symmetry and can 
be restored only in the continuum limit. With domain wall fermions, chiral symmetry 
can be effectively restored at finite lattice spacing, i.e. before taking the continuum 
limit, albeit at some increased computational cost. 

In chapter ^ after a brief overview of continuum and lattice QCD formulations 
intended primarily to fix the notation used in the rest of the work, we then discuss 
the domain wall fermion formulation and the related overlap fermion formulation. 
In particular, we develop a framework, based on the Gell-Mann, Oakes, and Renner 
(GMOR) relation, to discuss how domain wall fermions restore chiral symmetry on the 
lattice. Finally we discuss in general terms the physics underlying the deconfinement 
phase transition of QCD. 

In chapter |21 we discuss some specific technical issues regarding the efficient 
implementation of domain wall fermions for numerical simulation with dynamical 
fermions. We discuss the even-odd preconditioning of the fermion matrix. We also 
derive the molecular dynamics equations of motion for the fermion action for use in 
the hybrid Monte Carlo (HMC) algorithm. 

In chapter 121 we present the results of numerical simulations of Nf=2 QCD with 
light degenerate quarks. Because domain wall fermions are computationally (9(10)- 
0(100) times more expensive than standard lattice fermion formulations, we focus 
our efforts on studying the physics of the finite temperature transition region where 
the role of chiral symmetry should be of primary importance. Our initial studies 
on 8^ X 4 lattices located the transition region. Further studies on 16^ x 4 lattices 
searched for conclusive evidence of critical behavior in the transition region which 
might indicate the order of the transition. Finally, we make a quantitative estimate 
of the degree of chiral symmetry restoration in the broken phase using the GMOR 
framework introduced in chapter 
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We use the following conventions throughout this work. All lattice fields and 
bare parameters are to be dimensionless. To compare these dimensionless lattice 
numbers with their physical dimensionful analogues, we employ the dimensionful 
lattice spacing. For example, the lattice pion mass at a given lattice spacing a 
may be compared to the physical pion mass rni^^^'' by mf^^^'' — mT^j a. 
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Chapter 1 

QCD with domain wall fermions 

1.1 Continuum description 

QCD, or Quantum Chromodynamics, is the current best description of the strong 
nuclear force. The theory contains six types, or flavors, of fermions called quarks. 
Analogous to the theory of electrodynamics, where fermions called electrons carry 
electric charge, quarks carry one of three types of color charge. As electrons interact 
through an electromagnetic gauge field, these color-charged quarks interact through 
color-electric and color- magnetic gauge fields following a logical extension of Maxwell's 
field equations. The complex charge structure of QCD yields a far richer theory 
with quark confinement and gauge field self-interactions not found in the simpler 
electro dynamics . 

QCD is described schematically by the Lagrangian density 

^ -^quarks "l~ -^gluons ( • ) 

which yields the action when integrated over space-time. By treating the time coordi- 
nate as a complex parameter, we can equivalently define the Lagrangian either on the 
real time axis {Minkowski space-time) or the imaginary time axis {Euclidean space- 
time), as the theory defined on one can be analytically continued to the other. We 
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choose the Euchdean prescription because it allows the interpretation of the action 
as a Boltzmann weight in an equivalent statistical mechanics framework. 

The first term in the Lagrangian density includes a description of how quarks 
interact with gluons, the bosons that constitute gauge fields 

Nf 4 

>Cquarks = ^^'^j [if^ (^M " ^^^m) + ^j] ^Ij ■ (1-2) 
j = l tM=l 

The 7^ are Hermitian matrices which satisfy the anticommutation relations {7^, 7i/} = 
25/i,iy appropriate for Euclidean space-time. Nf is the number of quark flavors (six), 
the gauge fields are elements of the algebra of the Lie group SU(Ai'c) 

= E (1-3) 

and Nc is the number of color charges (three). The A° are the generators of the 
algebra of SU(A^c) and are normalized 

tr (A'^A'') = 25"^ (1.4) 

The second part of the Lagrangian includes a description of the self-interaction 
of the gluonic fields in the absence of quarks 

■^gluons = '^^,i/=l 2'^^ [-P'fiuF^u] ■ (1-5) 

F,u = Eali' = {d,A: - d^Al + gr'-^AlAl) f (1.6) 

where F^y is the Yang-Mills field strength tensor. Note that the summation convention 
is implied in the right hand side of ()1.6p . are the structure functions of SU(Ai'c) 
and satisfy the commutation relations 

[A", A^] = 2zr'"=A^ (1.7) 
The action is the Euclidean spacetime integral of the Lagrangian density 

S (A, q,q) = j (fx (^quarks + >^^gluons) • (1-8) 
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The partition function is given by the path integral over all possible field configura- 
tions 

Z = J [VA][Vq][Vq] e'^^^-^'^). (1.9) 

Physical quantities are computed as a weighted average of appropriately chosen ob- 
servables O over all field configurations, normalized by the partition function 

{0{A,q,q)) = [VA][Vq][Vq] 0{A,q,q) e'^^^'^'''). (1.10) 

As written, this partition function is formally divergent because it includes integration 
over an infinite number of gauge configurations which are related by gauge transfor- 
mations, thus of equal weight in the path integral. Restricting the path integral to 
sample only one element of each set of gauge-equivalent configurations will resolve 
this problem. For example, the method of Faddeev and Popov can be used 
Even still, there are further divergences which must be carefully regularized to ensure 
the partition function is well defined. 

In addition to invariance under local SU(A'"c) gauge transformations, when the 
mj=0, the action is also invariant under independent global transformations of the 
left and right chiral components of the fermion fields. We define the chiral matrix 
75 as the last remaining linearly independent matrix that anticommutes with the 7^ 
matrices: {75,7/x}=0. We then define the chiral projection operators 

p.-'-^ (1.11) 

and the right and left components are given by q^ = P+q and = P_g. Now the 
transformations explicitly are 



qn ^ e*"«Vg^, 



qR ^ e^'^qn, 



qL e'^'-qL. 
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Formally, these transformations generate the symmetry group SUl(A^/) ® SUr(A^/) ® 
Ui(l) C?> U/j(l). Furthermore, we stress that these are the symmetries of the classical 
action only. The full theory will only respect a subgroup of these symmetries due to 
quantum interactions. 

The Noether currents commonly associated with these symmetries are the non- 
singlet vector current 

f^{x)=q{x)^,^q{x), (1.12) 
the non-singlet axial vector current 

A" 

iM5(^) =^(^)7M75yg(a;), (1.13) 

the singlet vector current 

Jf^ix) = q{x)-ff,q{x), (1.14) 
and the singlet axial vector current 

J^^) = q{xh^i5q{x). (1.15) 

The conservation of the singlet axial vector current ()1.15|) is broken by the Adler- 
Bell-Jackiw (ABJ) anomaly. The conserved charge associated with the singlet vector 
current p.l4j) is the baryon number and is conserved even when the bare quark masses 
are non-zero. 

For massless quarks, the remaining non-singlet chiral symmetry may be sponta- 
neously broken by the ground state: SUiiNf) ® SUij(A^/) SUy(A^/). This implies 
that the non-singlet vector current ()1.12j) is still conserved. It remains conserved even 
for non-zero quark masses, provided they are degenerate, and the associated conserved 
charge is called isospin. The implied breaking of the non-singlet axial vector current 
is responsible for the appearance of the N"^ — 1 Nambu-Goldstone bosons. Though 
quark mass terms explicitly break this current, for small masses spontaneous and 
explicit symmetry breaking effects work together to produce light Nambu-Goldstone 
bosons whose mass is proportional to the square root of the quark mass. 
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Another interesting consequence of chiral symmetry is the role it plays in the 
spectral representation of the Dirac operator. From equation (jl.2j) . we define the 
continuum Dirac operator 



This anti-Hermitian operator has a continuous imaginary eigenvalue spectrum: Pqx = 
iXqx, where each eigenvalue occurs with the density p (A), normalized so J dX p (A) = 
1. Because p anticommutes with 75, then for each non-zero eigenvector qx we can 
generate a new eigenvector 'j^qx = q-x whose eigenvalue is the negative of the original 



Thus, non-zero eigenvalues always occur in pairs of opposite sign: (zA, 

A simple application of the spectral representation is to derive the spectral iden- 
tity that relates the chiral condensate to the pseudoscalar susceptibility in QCD with 
degenerate quarks. The chiral condensate in a volume V is 



where we have inserted an extra minus sign so that the symbol (gg) always represents 
a positive number for positive m. We have chosen the normalization such that (gg) ~ 
1/m as m ^ 00. p(A), an even function, is the quantum average of the eigenvalue 
density so it depends implicitly on the bare masses and couplings. 

The pseudoscalar, or pion, susceptibility is (no sum on a) 
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where we have imphcitly assumed that Nf > 2. We can anticommute P ^ and 75 

X.. = j^Tl'((m-P)-'(m+p)-'> 

^~,A^^g) . ,1.20) 



00 



(m — iX){m + iX) 

If we multiply the right hand side of (jl.l8p by the identity factor (m — iX) / (m — iX) 
and note that zAp(A) is an odd function, then (jl.lSp and ()1.20|) give the identity 

{qq)=mxn-. (1.21) 

Finally, we emphasize it is the chiral symmetry, manifested in p.l7|) . that is respon- 
sible for this identity. Not surprisingly, (jl.21|) can also be derived as the integral form 
of the flavor non-singlet axial Ward-Takahashi identity, another consequence of chiral 
symmetry. In section [T31 we will show that ()1.2H) . within the context of the spectral 
representation, provides a useful starting point for understanding chiral symmetry 
breaking for lattice fermions. 



1.2 Lattice actions for gauge fields 

Lattice gauge theories provide a regularization scheme for the divergences of the 
path integral by discretizing the Euclidean spacetime manifold. Not surprisingly, 
this discretization process destroys the continuous Lorentz symmetries, leaving only a 
hypercubic remnant. Other symmetries might also be broken if we make a poor choice 
when transcribing our continuum action into the lattice scheme. Wilson demonstrated 
j3] that the local gauge symmetry, at least, can always be exactly preserved on the 
lattice and other symmetries are restored in the limit of vanishing lattice spacing: 
a ^ 0. 

We can define the lattice partition function as the integral over all field configu- 
rations of the lattice action 

Z = [ [DU^] [Vq] [Vq] e-(^G+5^)_ (i_22) 
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Purely gluonic contributions to the action are contained in So- Fermionic contribu- 
tions are included in Sp and will be discussed in section IT751 

Starting from the A^{x) fields of (jl.3|) . we define the gauge links U^{x) as the 
integral of the path ordered exponential of the ^4^(2;) along the link connecting the 
site X to the neighboring site one lattice spacing away in the /i direction: x + fi. As 
the A^{x) are elements of the algebra of the gauge group SU(A^c), the gauge links are 
elements of the group itself. Additionally, the U^{x) are used to parallel transport 
quark fields from the site x + fito the site x. 

The Wilson action for SU(Ai'c) gauge theory in d dimensions is [3] 



When the plaquette is expanded in powers of the lattice spacing, the leading term 
appearing at Oi^a?) corresponds to the continuum gauge field Lagrangian density 
(ll.5|l . Terms in the expansion which appear at higher powers of a are called irrelevant 
because their infiuence vanishes in the continuum limit. Although the analogy with 
common statistical mechanics notation is intentional, here /3 is a bare parameter of 
the lattice theory which corresponds to the gauge coupling g in the continuum theory: 
(3 = 2Njg^ as a ^ 0. 

Wilson later proposed j3] introducing irrelevant operators to the gauge action to 
improve the scaling of Monte Carlo simulations. Here we consider adding only the 
rectangle operator 




(1.23) 



where /U, G [1, d\ and the plaquettes Up{x, /i, u) are 



Up{x, /i, u) = U^{x)U,{x + fi)Ul{x + v)Ul{x). 



(1.24) 




(1.25) 
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We require that the coefficients satisfy the constraint cq + 8ci = 1, as predicted in 
weak couphng lattice perturbation theory. The rectangles f//j(x, /x, z/) are 

Ur{x, /i, u) = U^{x)U^{x + fi)Uy{x + 2fi)Ul{x + fi + v)Ul{x + v)Ul{x). (1.26) 

In chapter ini we will compare numerical simulations that include a different choice of 
gauge action looking for behavior which appears closer to the continuum limit. 

1.3 Domain wall fermions on lattice boundaries 

The standard approach to discretizing the continuum quark Lagrangian density p.2|) 
is to define the quark fields on lattice sites and replace derivatives with finite dif- 
ferences. For example, we can define the gauge covariant lattice central difference 
operator corresponding to the covariant derivative in the continuum 

V(/i)^,^, = ^ [Ui,{x)5^+(,^^. - UjXx - /i)4-A,^'] • (1-27) 

As we shall see below, this naive prescription leads to the well known problem of 
species doubling: on a c? dimensional finite lattice the discretized action describes 2*^ 
identical species, one for each corner of the Brillouin zone. Nielsen and Ninomiya 
13 El showed when the lattice fermion action is local, Hermitian and translationally 
invariant then either the fermions are doubled or chiral symmetry of the continuum 
action is broken on the lattice. 

Wilson proposed giving all but one of the doubled fermions a mass proportional 
to the lattice cutoff by adding an irrelevant term to the fermionic action. Using the 
lattice central difference operator corresponding to the Laplacian in the continuum 

Mf^)x,x' = [U^,{x)6^+f,y - 26^y + Ul{x - , (1.28) 

the lattice action for Wilson fermions is 
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We suppress spin and flavor indices and the summation convention implies summing 
x,x' over all sites on the lattice. The bare fermion fields ip and if), the bare mass 
mo and coefficient of the irrelevant term r are all dimensionless quantities whereas 
V(/i) and A(yu) are dimensionful according to the powers of the lattice spacing in 
their definitions. Rearranging terms, we put the action into the standard form of a 
fermion matrix bilinear 

Sr=i^.MZh.' (1.30) 

where M^^} is the Wilson-Dirac fermion matrix 

d ^ 

(1.31) 

For numerical simulation, r is usually set to one to take advantage of the simplified 
spin structure of the spin projection matrices: P±^ = (1 ± 7^) /2. 

The free Wilson fermion propagator is obtained by inverting the fermion matrix 
with all Ufj_{x) set to the identity. We block-diagonalize the fermion matrix by Fourier 
transform to give 

d 

M^^) = mo - [27^ sinp^ + r (cosp^ - 1)] , (1.32) 

where the are dimensionless momenta that reside on the reciprocal lattice and 
range over — tt < < tt. In the massless limit r — > corresponding to the naive case 
above, we see that the propagator has a massless pole when each = or tt. These 
are the corners of the Brillouin zone mentioned above and each pole is associated with 
one doubler species. When < r < 1 only the fermion species with p^ = remains 
massless and the others acquire a mass on the order of the lattice cutoff. However, 
even though we have removed the doublers with the irrelevant Wilson term, chiral 
symmetry is not restored because the Wilson term commutes with 75. 

The domain wall fermion formulation, originally proposed by Kaplan JH], solves 
the doubling problem while conserving the SU L{Nf) (g) SU R{Nf) chiral symmetry 
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of massless quarks on the lattice through the introduction of an infinite internal 
dimension which we label with coordinate s. To demonstrate the construction, we 
start with an odd d+1 dimensional action of free lattice Wilson fermions, see ()1.29|) . 
where the bare mass parameter mo is a function of one direction, the s-direction. 
Although mo(s) can be any monotonically increasing function that crosses zero, for 
simplicity we discuss only step function mass terms 

{m+, s > 
(1.33) 
m_, s < 

Despite the presence of a mass term and the lack of chiral symmetry in odd dimen- 
sions, there is a normalizable, chiral basis of states that solve the equation 



p 



^27^ sin 



vi/+ (1.34) 



where M*^^^ is Kaplan's Dirac matrix and is a d dimensional momentum. The 
existence of such a basis indicates the spectrum contains states that propagate as 
massless fermions in the even d dimensions. Kaplan's solution 

^+ = (±l)^e*P-^e-^°l'l, sinh/io = mo (1.35) 

indicates the massless modes are exponentially localized on the domain wall, where 
the mass changes sign. Furthermore, when the extent Ls along the s-direction is 
infinite, only one chirality of the massless modes is normalizable. So, even when the 
lattice spacing is finite, the action describes a single massless chiral fermion flavor at 
scales well below the lattice cutoff in the ^ oo limit. 

On a lattice with a finite Ls, the solutions of the opposite chirality become nor- 
malizable and localized on the boundaries. In essence, the boundaries become anti- 
domain walls and at low energies the effective theory is naturally vector-like. As 
the separation between the opposing walls is increased by increasing Lg, terms that 
break chiral symmetry will eventually become negligible due to the exponentially sup- 
pressed overlap between modes of opposite chirahty on opposite walls. So, domain 
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wall fermions at large Ls are ideally suited for simulating the massless quarks of a 
vector-like theory like QCD. 

If we reinterpret the s-direction as an internal flavor space, then gauge fields are 
coupled to the fermions in the natural way [HI HDl HI]- We write this as the set of 
constraints on the gauge field 



Ui,<d{x, s) = U^{x); Ud+i{x, s) = 1. 



1.36) 



From this perspective, domain wall fermions are identical to heavy Wilson fermions 
in d dimensions that carry an additional flavor index s and derivative terms along 
this s direction are now interpreted as flavor mixings in a sophisticated mass matrix. 

Shamir proposed J2] a variant of domain wall fermions where free boundaries 
at both ends of the extra direction serve as the domain wall/anti-wall. This variant 
naturally arises when considering the step function mass term ()1.33p in the limit 
m_ ^ — oo, provided we use some care in handling the boundaries. From a practical 
view, this reduces by half the number of lattice sites necessary to separate the walls 
a fixed distance. The action is a fermion bilinear 



(-.(dwf) _ TjT j^Adwf) J 



;i.37) 



Starting with the Wilson fermion matrix p.31|) we build the flavor mixing matrix 
from first and second order difference operators and the chiral matrix = 7^+1 of 
the effective even d dimensional theory 



2 



:i.38) 



The hat notation on the difference operators denotes that the operators are ungauged 
and dimensionless because both act on a flavor space: 



^s,s' = ^s+l,s' ~ 25s, s' + ^s-l,s 



;i.39) 
;i.4o) 
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Following the literature, we collect terms in the action that are diagonal in the 
s-direction into an operator that acts parallel to the d dimensional spacetime 

^2l' = ^L'^.,.' + 5.,.'^t'- (1-41) 

The parallel piece is simply a Wilson matrix with slightly different mass term, d\ ^, = 
^^i* + rb^x,x'-, which we write out explicitly for completeness 

dI,x' = i^d + rs - mo) 6^,^^' 
d ^ 

- 5Z 9 [(^ - ^A^) U^{x)5,+i.,x' + (r + 7m) Uli"" - A)^x-A,x'] • (1-42) 
M=i ^ 

The second term is the flavor mixing matrix (with flavor index s) 

Dts'=\ -(^)W , 0<.<L,-1 . (1.43) 

We have included a new parameter rrif that plays the role of generalized boundary 
conditions in the s-direction. For rrif = the boundaries are free and for = ±1 
and rs = 1 the boundaries are (anti)periodic. We shall also see that ruf controls the 
bare mass for the effective light modes. 

Until now, we have included a parameter analogous to the Wilson r because it 
is possible to find massless modes in the low energy theory for < rs < 1 8J. While 
r plays the physical role of setting the lattice mass scale of the 2^^ — 1 doubler states 
for Wilson fermions, has no such physical interpretation. Furthermore, for rs 7^ 1 
the massless modes are no longer eigenvectors of Fs. So, for the remainder of this 
work we will assume rs = 1. 

As in fll.32j) . we can diagonalize Z?" in the free case by Fourier transform 

d 

5| = b{p) - i-ff, sinpf, (1.44) 

M=i 

where b{p) = 1 — mo — X]m=i ^ (cosp^ ~ !)• As in our previous discussion of Kaplan's 
formulation, we want to search for a normalizable basis of states upon which the full 
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domain fermion matrix acts like the chirally symmetric Dirac matrix in d dimensions 

d 

^p^J'^P.^' = - 5Z '^/^ sinp^^p,,. (1.45) 

Since we chose = 1, the \1/ should also be eigenvectors of T^. = 4>p^sU^ where 
are constant eigenvectors of T^. 

We can reformulate the problem of identifying this basis into solving the zero 
eigenvalue problem of a different matrix 

d 

Mp;s,s' ^ M^tJ + E sinp^^,,,, = b{p)Ss,s' + Dj^,,. (1.46) 

Note that Ai essentially contains all the terms in M*^'^™^-' that can break chiral sym- 
metry in the d dimensional spacetime. So, the basis of chiral states corresponds to 
the null space of Ai 

Mp.,s,s'^p,s' = [b{p)S,,s' + Dj^,,] (P+</)+,, + P_0-,,) = (1.47) 

where P± = (1 ± T^) /2 are chiral projection matrices analogous to ()1.11|) . For mj=0, 
the eigenvectors which solve this equation are 

If we can choose mo to satisfy the condition 

d 

-1 < l-mo-^(cosp^-l) < 1 (1.49) 
^l=l 

then is localized at the boundary s = and, as — ^ oo, it is the only normalizable 
solution at any finite s. For the zero momentum state, < mo < 2 will satisfy the 
condition. Similarly, for 2r < mo < 2 (1 + r), d zero momentum doubler states will 
be bound to the domain wall. It is straightforward to derive the other conditions for 
the remaining doubler states. If we allow Ls to be finite, becomes normalizable 
and localized at s = — 1. 
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We remind the reader that the states \l/p^s are not eigenvectors of the free do- 
main wall fermion Dirac matrix. Since that Dirac matrix is non-Hermitian, it is 
essential to calculate the free propagator from the Hermitian second order operator: 
]\^t(dwf) j\^(dwf)^ Por non-zero mj, Vranas ^3 E] analytically computed the wave 
functions of the nearly massless modes and found the effective mass (for r = 1) is 



The Ls-dependent contribution is exactly the overlap between the wave functions of 
the right and left modes on opposite walls. 

In the interacting case, the precise shape of the wave functions of the nearly 
massless modes depends upon each gauge configuration. So it is convenient to define 
the effective light fermion fields through the projection 



For large Lg, these fields have a large overlap with the nearly massless modes and a 
vanishing overlap with the other heavy modes described above. 

Besides the nearly massless modes, the action also describes 2*^ — 1 fermion dou- 
blers that propagate along the wall as massive fermions with the mass of 0{r/a). Fur- 
thermore, naive counting suggests there are an additional 2'^ {L^ — 1) massive states 
in the spectrum. These states will contribute to a bulk infinity as — oo. In section 
11.41 we will identify this divergence and introduce terms into the action to cancel it 
and preserve the — > oo limit. 



1.4 The limit of infinite domain wall separation 






Using the rules for integration of Grassmann variables, we can explicitly perform 
the integration over the fermionic degrees of freedom in the lattice partition function 
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()1.22|) . The resulting expression now depends explicitly on only the gauge degrees of 
freedom 

Z = j [Pf/^] (det M) e-^« = j [VUf,]e-^^'' (1.52) 

where M is the lattice Dirac fermion matrix which appears in the fermion action, e.g. 
(I1.31|) or (jl.38p . and Scs = Sq — TrlogM. Since the fermionic weight in the path 
integral is given by determinant, we can consider replacing M with a different matrix 
M' provided det M = det M'. As we shall see, this freedom allows us to uncover 
important details about the physics of the lattice fermions. Following the approach 
of Neuberger ^3] and Kikukawa and Yamada ^Hj, we will construct a product of 
matrices whose determinant is the same as the domain wall fermion Dirac matrix in 
section [T731 yet whose Lg oo limit is manifestly apparent. 

First, we choose a Weyl basis for the Dirac matrices 

l,=i,-4 = 1 ° (1-53) 
where each block is a 2^^^ x 2^^^ matrix. We define the q x q matrices 

d 

2 



B = {rd + l-mo)Sx,x'~7:^[Uf,ix)6x+f,,x' + Ul{x-fi)6x-f,,x'], (1-54) 



d 



C = ^ o'm \-Ufi{x)6x+ii,x' + Ujjx - fi)Sx-f,,x'] ■ (1-55) 



2 



where q = 2^~^NfNc^ and is the number of lattice sites in the d dimensional 
spacetime. This block notation allows us to write the matrix Z^" of section [T31 see 
()1.42|1 . in the compact form of a 2g x 2q block matrix 




(1.56) 

Continuing with the same strategy, the full domain wall fermion Dirac matrix in 
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()1.41|) can be written as the block matrix 



]\/[(<iwf) 



B -C 
Ct B 
-1 





-1 



m/ 




v 





m/ 



-1 






-1 

B -C 
Ct B 



1.57) 



It is a potential advantage of this block form that we can permute blocks of 
rows and columns, leaving the determinant unchanged up to a sign. If q is even, 
which is always true if > 4, then the sign is always positive. We choose the set of 
permutations to make M^'^'"-^^ as block lower triangular as possible, leaving a single 
non-zero block in the upper right corner. This is done by swapping pairs of block 
columns and then moving the first block-column on the left all the way to the right. 
The new matrix M''^'^"'-^) in terms of 2q x 2q blocks is 



/3i 



\ 



\ 



:i.58) 



where the and l3s are block lower triangular 

B 

Ct -1 

-1 -c 

B 




;i.59) 



;i.6o) 



We use the technique suggested by Neuberger to simplify the determinant of 
the nearly lower triangular M'^'^^f\ First, we factor (|1.58|) into the product of a lower 



21 



triangular and an upper triangular matrix 



a 



V 



a 



( 



;i.6i) 



where the Vg are given by 



-VTLf 



1.62) 



1 
J 

As a is block lower triangular, it is straightforward to compute it's inverse in block 
form. Now the determinant is 



-nif 
1 



-1)"^= (detS)^^ det 



T 



T- 



det M('^'"^) = (det ay-'^ det ai, det (1 - vlJ (1.63) 

1 y y -m/ 

where it simplifies what is to follow if we define the 2q x 2q matrix T and its inverse 

B-^ -B-^C^ 
-CB-^ CB-^C^ + B 

C^B-^C + B C^B-^ 
B-^C B-^ 

At this point we have already made substantial progress since the dimensionality of 
the matrices on the right hand side of ()1.63p no longer depend on L^. 

However, to fully expose the Lg ^ oo limit, we should find a way to raise to 
an arbitrarily large power. Let 1 be the 2q x 2q identity matrix and be the 2q x 2q 
block matrix 

^ 1 

15 ^ I I . (1.66) 

-1 



;i.64) 



;i.65) 
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Using these matrices and the identities 




[(i_^^)i + (i + ^^)r5] 



-[{l-mf)l-{l + mf)V,] 



;i.67) 



1.68) 



we can bring det M^'^'^^^ into the form 



det M^'^"'^) 



-1)"^^ (det Bf' det (l + T'^') det T. 



X det - 
2 



'1 + mf)l — (1 — mf)T^ tanh | — — logT 



1 + m/)l + (1 -m/)r5tanh(^i/) 



X det 



1.69) 



where in the second equation we have used detFs = (—1)'^ and defined the 2q x 2q 
matrix H = — log T. 

Finally, we see the Ls ^ oo limit of p.69|) is divergent from the factors (det B)^" 
and det (l + e^"^^. This should not come as a surprise if we recall that domain wall 
fermions can be interpreted as a theory of Lg heavy Wilson fermions that mix to 
create low energy chirally symmetric massless modes. We should remove this bulk 
divergence if our partition function is to correctly model a single light fermion fiavor. 

A simple method to remove the contribution of the heavy modes was proposed 
by Vranas E]- For rrif = 1 the last determinant factor in (jl.69j) is unity. So, if 
we divide the partition function for arbitrary mj by the partition function with fixed 
mf=l, we will cancel the bulk divergence and render finite the — > oo limit. The 
technique used to place the determinant into the denominator is to add Pauli-Villars 
regulator fields^ into the action with the same Dirac matrix as the fermion fields except 

with fixed mj=l. In the path integral over the Pauli-Villars fields, the integration 

^Pauli-Villars fields, also called pseudofermions, commute like bosons yet have the same matrix 
spin structure as fermions. 
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is now Gaussian rather than Grassmannian, so the determinant of the Pauh-Villars 
matrix will appear in the denominator. Explicitly, the proposed partition function is 



^(dwf) 



[VU] [VW] [V^] 



, det M(mf) ^ 



[V<^Pv] e' 



[VU] det 



detM(m/ = i; 

1 + mj)U + (1 - mf)T^ tanh(^i/) 



-Sg 



;i.70) 



where 



(1.71) 



is the Pauli-Villars action defined as in ()1.37|) . As we shall see in chapter |21 it is 
sometimes convenient to also rewrite the fermionic path integral in terms of pseud- 
ofermions as well 



^(dwf) _ / [p^J ^p^^j ^p^^J 



V^ 



t 

PV 



[V<!)pv] e 



1.72) 



where the fermionic action is now written as 



(1.73) 



To take the Lg ^ oo limit of p.70|) . we recall that the hyperbolic tangent function 
converges to the sign function 



lim tanh(nx) 



eix 



1, X > 
-1, X < 

Then the fermion determinant for QCD with exact chiral symmetry is 



1.74) 



det M(m 



f) 



detM(m/ = l) 



det 



l + m;)]l + (l-m;)r56(if) 



1.75) 



So, we have achieved our goal of finding the contribution to the partition function of 



domain wall fermions in the Lg oo limit. 



As an aside, it is possible to construct a transfer matrix formulation of domain 
wall fermions analogous to the Hamiltonian formulation of lattice gauge theory. The 
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crucial difference is that the domain wall Hamiltonian corresponds to evolution along 
the s-direction rather than in time. Due to the trivial gauge links in this direction, 
it is possible to construct a Hamiltonian that is s-independent. The matrix T in- 
troduced in ()1.64|) and its logarithm H are the s-independent transfer matrix and 
Hamiltonian of this formulation. Finally, we note the domain wall fermion transfer 
matrix formulation described here is just a specific case of the overlap formulation 



In the preceding section, we have described how domain wall fermions restore chiral 
symmetry in the limit of infinite domain wall separation. However, for numerical sim- 
ulations we must necessarily work at finite domain separation. For physical quantities 
defined at large distances, we expect that the effect of the chiral symmetry breaking 
induced by finite Lg is to shift the bare mass by some small additive piece m-^esi 
the residual mass. By examining the non-singlet axial Ward-Takahashi identity, we 
hope to identify and measure the residual mass in the confined phase of QCD. 

The domain wall fermion action ()1.H7|) with degenerate quark masses is invariant 
under a global SUy(A^/) symmetry. Following Furman and Shamir J7j, we can define 
the corresponding d+1 dimensional conserved non-singlet vector current as 



imiiHiiin]. 



1.5 Effects of finite domain wall separation 




< s < - 



s = 



1 





We can see immediately this current is conserved as its divergence vanishes 



d+l 
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d 



However, as the low energy effective theory is d dimensional QCD, we can also define 
a conserved d dimensional vector current. Towards this end, we move the terms in 
()1.78|) for /X = d + 1 to the right hand side to define the continuity equation 

-JS+i{^, 0) - - 1), s = 

[ rrifj^^.ix, - 1) + J^^,{x, - 2), s = - 1 

(1.79) 

where d^J^{x,s) = J^l^{x,s) — J^{x — fl,s) for yU = l,---,d. If we define the d 
dimensional non-singlet vector current as 

Ls-l 

v;(x) = 5^ (1.80) 

s=0 

then we can see from ()1.79j) that it is also divergenceless: X])!=i '^AtH"^^) ~ 
isospin is a symmetry of the low energy effective action of domain wall fermions. 

To proceed as above to construct a d dimensional axial vector current is not so 
straightforward as it cannot be uniquely defined. Clearly, opposite chiral rotations 
should be applied to the states on opposing walls. Following Furman and Shamir 
fri\ we perform opposite chiral transformations on each half of the fermions in the 
s-direction. Analogous to ()1.80|) . the corresponding d dimensional non-singlet axial 
vector current is 

where our notation assumes, for simplicity, that Lg is even. Unlike the vector current 
in p.78|) . the divergence of the axial vector current has two contributions that do not 
trivially vanish 

d 

^9^^;;(x) = 2m;J,'Vi(x,L, - 1) + 2J,'Vi(x,L,/2). (1.82) 

When expressed in terms of the fields qx,qx of f)l-51|) . Jl^^i{x, L^ — 1) takes the familiar 
form of the pseudoscalar density 

J-;+i(a;,L,-l) = g,75yg.. (1.83) 
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Thus, when ruf ^ the first term reproduces the correct contribution expected from 
a naive fermion mass term. Although the second contribution, J'^_^^{x, Ls/2), does 
not vanish for finite Lg, we present numerical results in section that indicate it is 
exponentially suppressed. 

To better understand the consequences of working at finite Lg, we will follow 
the role of J'^j^i{x, Ls/2) as we derive the Gell-Mann, Oakes, and Renner (GMOR) 
relation PU] in the domain wall framework [^|221- We start with the axial Ward- 
Takahashi identity 

^ I A" \ / A" A" 

+2^J,%(x,L,/2) go75ygo| 
-Sxfi (go9o) (1-84) 

where the last term comes from the derivative of the Heavyside function used for 
time ordering and the zero subscript of (^qQ'o) indicates the contribution to the chiral 
condensate from a single point, analogous to setting x=0 in p.l8|) rather than aver- 
aging over the whole spacetime volume. The left hand side of ()1.84|1 will vanish when 
summed over the whole lattice, leaving us with the GMOR relation 

{%qo) = mfXn,o + AJ5 (1.85) 

where Xn is the pion (or pseudoscalar) susceptibility 

2 / A" A'' 



Xtt,0 



AN, 

and AjTs is the mid-plane contribution due to finite L 



Yl \ 1xl5^(lx go75ygo ) (1-86) 



X 



As before, the zero subscript Xvr.o in p.86|) means this quantity is computed only at 
a single lattice point, analogous to setting a;=0 in ()1.19j) . 
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In the chirally broken phase, both susceptibihties should be dominated in the 
massless hmit by the hght Goldstone degrees of freedom locahzed on the walls. As 
J^^-^{x^Ls/2) is located midway between the two boundaries, its overlap with the 
boundary states should be exponentially suppressed. So, assuming an infinite lattice 
volume, we will model Xtt,o 

Xn,o = — + ao + 0{mf) (1.88) 

and AjTs will also have a pole when mf = — mres although we expect that the q will 
be exponentially suppressed in Lg 

AJ5 = — + Co + O{mf). (1.89) 

rrif + mres 

Thus, — > as ^ —mres- We also note that the parameters used here are 
assumed only to be independent of and can vary with ttlq and Lg but we imagine 
performing the following derivation while holding mo and Lg fixed. 

We write the GMOR relation as 

= m, + ^ (1.90) 

Xtt,0 Xtt,0 

and substitute our parameterizations ()1.88|) and ()1.89p in the right hand side 

(gpgo) _ ^ I c,i + Co(my + m,es) 

As we tune mj to the pion pole, i.e. rrif —mres, the left hand side of 1)1.911) 
will vanish. On the right hand side, this requires mres = c_i/a_i which allows us to 
eliminate c_i as a free parameter in ()1.91|) . 

In the continuum identity the slope of {%qo) /Xn,o versus the bare quark mass is 
always one. However, for domain wall fermions as rrif ^ —fnj-cs at fixed mo and Lg 



d ( (go9o) 



^^co^aomres ^^_g2) 



mf=—mrcs,mo,Ls ^ 

SO the slope is never one except when co = ao'^res- To clarify this situation, let us 
eliminate the parameter co in favor of bo = cq — aorrij-es, the parameter which, when 
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non-zero, indicates the slope of (go?o) /Xtt.o is not unity as rrif —>■ —rrires- In this new 
parameterization, we see that ()1.89j) becomes 

AJs = mresXTT.o + bo (1.93) 
and the GMOR relation becomes implicitly 

{QoQo) = il^f + ^res) Xn,0 + (1-94) 

and explicitly 



ho 



a_i + ao (m/ + nires) 



;i.95) 



Finally, we recall that the chiral condensate receives contributions from chiral 
symmetry breaking not just at large distances but from all energy scales up to the 
lattice cutoff. We should expect that chiral symmetry violations of the GMOR rela- 
tion due to lattice artifacts will appear in a form that cannot be interpreted as an 
additional residual mass in the low energy effective Lagrangian. We conjecture that 
the non-zero 69 is precisely the cumulative contribution due to the chiral symmetry 
breaking of the heavy modes at the lattice cutoff. As such, we could absorb the 
heavy contribution into a new definition of the chiral condensate, which we call the 
"subtracted" chiral condensate 

(%9o)sub = (^o9o) - &o- (1-96) 

One potential advantage of including this contribution in the chiral condensate is that 
it guarantees the slope of the lattice spectral identity will be unity when the spectrum 
becomes chiral as — —m^cs at fixed mo and Lg 

mo, is 



dnif \ Xtt.o 



In chapter El we will use the approach of this section to estimate the residual 
mass for our simulations near the critical temperature of QCD. For a given domain 
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wall separation, we would like to establish whether the dominant contribution to the 
effective quark mass comes from the bare input mass, rrif, or the residual mixing of 
the chiral surface states, rrij-es- Ideally, to ensure good control over the chiral limit, we 
would prefer to choose Lg large enough so mres rrif. Only then can we be confident 
that any observed critical behavior, or lack thereof, is a consequence of the chiral 
symmetry of the QCD Lagrangian. 

1.6 The QCD finite temperature phase transition 

An outstanding success of lattice QCD has been the demonstration of confinement 
and spontaneously broken chiral symmetry at low temperatures and the existence 
of a phase transition above some critical temperature, T^. A current goal of lattice 
theorists is to calculate the physical value for Tc, as well as other physically observable 
consequences of the transition, for comparison with anticipated experimental results 
from the Relativistic Heavy Ion Collider (RHIC) and other experiments designed to 
study the thermodynamics of QCD. 

In either the quenched limit with infinitely massive quarks or the chiral limit 
with massless quarks, the finite temperature phase transition is characterized by the 
behavior of an observable which is not invariant under some symmetry of the action. 
In the symmetric phase, the observable has a vanishing expectation value because 
the ground state of the system respects the symmetry of the action. In the broken 
phase, the observable has, in general, a non-vanishing expectation value because the 
ground state of the system no longer respects the symmetry. The observable is an 
order parameter for the transition. 

For quenched QCD, the gauge actions (|1.23|) and (|1.25|) are globally symmetric 
under a multiplication of all time-like gauge links originating from a given space-like 
hyperplane by the same element of the center of the gauge group. For SU(A^c), the 
center is Z^^. However, a closed loop of link variables with non-trivial winding around 
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the temporal direction of the lattice, commonly called a Wilson (or Polyakov) line, 



is not invariant under the transformations described above since it involves only one 
time-like link originating from each space-like hyperplane. The Wilson line observable 
measures the free energy of a isolated static fermion which couples to the gauge fields 
in representation R of the gauge group 1211 1211 12^1 • In the low temperature phase 
of quenched QCD, the free energy should be infinite and the expectation value of the 
Wilson line zero. Hence, the Wilson line is the order parameter for the quenched 
QCD phase transition and the symmetric phase is at low temperature. 

Lattice QCD studies have firmly established that the quenched transition is a first 
order phase transition. This means that a first derivative of the log of the partition 
function is discontinuous in the transition region. This discontinuity in the energy 
density is called the latent heat and is the measure of the amount of energy which 
must be put into the symmetric state to drive it into the broken state. 

For QCD with dynamical fermions, the Wilson line is no longer an order pa- 
rameter because the global Zjv^ symmetry of the gauge sector is not a symmetry 
of the fermion action. The physical interpretation is that once dynamical fermion- 
antifermion pairs are present in the vacuum, the static fermion can be screened by 
the vacuum, so the free energy must be finite. In the low temperature phase, since 
one doesn't expect to see isolated fermions the Wilson line expectation value should 
still be small. Thus, the Wilson line should undergo rapid evolution through the 
transition and can serve as an approximate order parameter. 

In the limit that the Nf > 2 dynamical fermions of QCD become massless, the 
appearance of the SUL(A^/)(8>SUij(A^/) global chiral symmetry discussed in sect ions ITTT] 
and 1 1.51 suggests the possibility of a new order parameter for the transition. The chiral 
condensate ()1.18|) is not invariant under chiral transformations. So, if the ground state 
of the system is invariant under chiral transformations then the expectation value for 



1 



■L4-1 



dim(i?) 
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the chiral condensate must vanish. However, if the chiral symmetry of the action 
is spontaneously broken then the chiral condensate will, in general, have a non-zero 
expectation value. Contrary to the quenched theory, the symmetric phase for the 
massless Nf > 2 theory occurs at high temperature. 

The order of the transition depends strongly on the number of massless quark 
flavors. For Nf > 3, renormalization group analysis indicates there are no stable 
infrared fixed points, a necessary condition for the transition to be second order, i. e. 
the second derivative of the log of the partition function is discontinuous in the critical 
region [2Zl I2H1 • Numerical lattice QCD simulations confirm this conclusion as first 
order transitions have been observed in the massless limit of QCD with A^j=3,4. For 
Nf=2, the order of the transition may be determined by the fate of the anomalously 
broken Ua(1) symmetry in the high temperature phase. If the Ua(1) symmetry 
remains broken while chiral symmetry is restored, then Nf=2 QCD is expected to be 
in the 0(4) universality class, and models in that class have second order phase 
transitions [211 • However, if U^(l) symmetry is effectively restored in the finite 
temperature phase transition then the universality class changes to 0(2) ® 0(4) and, 
as other models in that class are known to have first order transitions, it is possible 
that the finite temperature transition for Nf=2 QCD will be first order j2Hl EDI El- 
Earlier attempts to address this question using staggered fermions [121 jHll IMl 
ESI EHl did not produce conclusive results. While staggered fermions have an exact 
Uv(l) ® Ua'(1) chiral symmetry, the U^/(l) is actually an unbroken subgroup of the 
flavor non-singlet chiral symmetry and not the 11^(1) of the continuum. At finite 
lattice spacing, the flavor singlet Ua(1) is actually broken at the classical level and 
it becomes quite difficult to disentangle the signal of anomalous symmetry breaking 
from the explicit breaking due to lattice artifacts. 

As Ls — > oo, Uyi(l) is a symmetry of the classical massless domain wall action. 
The domain wall Dirac operator can have exact zero modes [IP] and a well-defined 
integer index, even for gauge fields at strong coupling. We demonstrated this numer- 
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ically in earlier studies [HZl IHH IMj by observing the effects of fermionic zero modes 
produced by smooth topological gauge configurations for small bare quark masses and 
for Ls > 10. Further, these effects persisted even after the addition of a moderate 
amount of Gaussian noise to the classical gauge fields. From this we conclude that 
domain wall fermions are an ideal formulation to study the Ua(1) problem. In prelim- 
inary studies ^01 1^ 1321 E31 E] ; we have strong evidence that Ua(1) remains broken 
in the high temperature phase, but by such a small amount that it may be possible 
that the breaking is negligible so that the symmetry is effectively restored. It is a 
primary goal of this work to address whether numerical lattice simulations of Nf=2 
QCD with domain wall fermions give evidence of a first or second order transition in 
the massless limit. 

In the real world, the fermions of QCD are neither massless nor infinitely massive, 
so both the Wilson line and the chiral condensate serve only as approximate order 
parameters. As such, we expect that the values of these observables will undergo 
a rapid change as we adjust the bare parameters of the theory to move through 
the transition region. However, the width of the transition region is an important 
phenomenological number and can be directly computed on the lattice. 

In sections 13.21 and 13.31 we will measure the Wilson line and chiral condensate in 
lattice simulations of Nf=2 QCD with two light degenerate flavors of domain wall 
fermions. We can locate the region of the finite temperature transition by searching 
for simultaneous rapid changes in both approximate order parameters as we vary 
the temperature. Once we have located the critical temperature, we will study the 
behavior of the chiral condensate for small quark masses in the hope that we will see 
some indication of the order of the phase transition in the chiral limit. 
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Chapter 2 

Numerical implementation of QCD 
with domain wall fermions 

2.1 General considerations 



A naive approach to numerical simulation of lattice QCD follows from our discus- 
sion in sections ll.ll and lL2| particularly equations ()1.10|) and ()1.22|1 . regarding the 
computation of expectation values. To evaluate the lattice partition function on a 
computer, we approximate the continuous multidimensional path integral over an in- 
finite lattice volume as a sum over finite volume lattice configurations with a finite 
precision representation of real numbers 



Each U is one possible lattice configuration and {U} is the set of all representable 
configurations. After fixing the lattice size and precision, only a countable number 
of configurations can be represented in the memory of a computer. So, to evaluate 
the expectation value in (j2.1|) . it is possible to program a computer to generate, one 
by one, each configuration U that can be represented, compute the observable 0{U) 






{u} 



{u} 
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and the Boltzmann weight e ^^^^ on each configuration and do the weighted average 
indicated. 

To see this approach as completely intractable, consider a simple Ising model: on 
each lattice site is a single degree of freedom which can take on one of two possible 
values, which we can represent as a single bit in the memory of our computer. If V is 
the total number of sites in a given lattice volume, then the number of representable 
configurations that our program must evaluate is 2^. The time required to compute 
any expectation value grows exponentially in the size of the configuration. For QCD, 
the only difference is the counting of representable degrees of freedom per lattice 
site once the number of bits of memory per lattice site is fixed. So, this method is 
intractable for all but the smallest lattice volumes. 

We recall from the statistical mechanics of systems with many degrees of free- 
dom, the Boltzmann distribution in configuration space is strongly peaked around the 
minimum of the action. This implies that the naive method above is incredibly ineffi- 
cient since most of the configurations generated will have vanishingly small weight for 
expectation values. To surmount this problem, we use importance sampling, which 
creates a sample ensemble of configurations drawn from the set of all representable 
configurations with a probability given by the Boltzmann weight. Configurations 
with large weights will appear more frequently in our sample. The number of con- 
figurations in the sample should be large enough that expectation values computed 
on the sample agree, within a given precision, with ones computed on the full set. A 
computer program could generate configurations at random and compute the weight 
and expectation values. However, this approach is no better than the naive method 
above since most configurations will be very unlikely to appear in a sample ensemble. 
So, an importance sampling method should provide us with an algorithm which only 
generates configurations with a good chance to be added to the sample. 

A classic example of an importance sampling method is the Metropolis algorithm. 
The algorithm has two parts: the update process and the acceptance test. The initial 
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condition for the algorithm is to start at any point in the configuration space, let's 
call it Ui. The algorithm then starts by computing the Boltzmann weight, Wi, of the 
initial configuration and adding it to the sample ensemble. Next, we use the update 
process to generate a second configuration with known probability P„ {U2 ^ Ui) and 
compute its weight, W2- Then we use the acceptance test to decide if the second 
configuration should be accepted into the sample with probability Pa{U2 Ui): 



The total probability that the second configuration is accepted given the first config- 
uration is P {U2 <— Ui) — Pa{U2 <— Ui) Pu{U2 <— Ui). We then iterate the algorithm 
by applying the update process and acceptance test to the last configuration accepted 
until the sample ensemble is sufficiently large. 

Two important issues must be addressed when designing update processes for the 
Metropolis algorithm. The first is whether the algorithm can eventually construct 
a sample ensemble and the second is how efficiently the algorithm can construct 
the ensemble. The first issue can be resolved if we establish two properties. The 
property called detailed balance, which guarantees that the Boltzmann distribution of 
our sample ensemble will converge to the distribution of the full set, requires 



Prom the preceding paragraph, it is straightforward to verify that the Metropolis 
algorithm satisfies detailed balance. The second property, called ergodicity, requires 
that the update process of the Metropolis algorithm can reach any point in the con- 
figuration space with finite probability in a finite number of steps. This property 
ensures that no important regions of configuration space are left unsampled provided 
the algorithm is run long enough. However, it is typically very difficult to determine 
analytically the efficiency of a given algorithm to construct an ergodic ensemble. In 
such cases, one attempts to address this issue through numerical studies. 




(2.2) 




(2.3) 
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For pure gauge theories on the lattice, a typical update process is to choose 
a single gauge link, propose a change to that link and compute the difference in 
the action, AS, between the two configurations. Since the gauge action is local, 
AS is inexpensive to compute. However, to maintain efficiency, some care must be 
taken in the update process to make sure that the proposed change to the single 
link is still sufficiently small to allow a reasonable chance of passing the acceptance 
test. Performing a Metropolis step at each link in the lattice once in succession is 
commonly called a sweep. One potential disadvantage of this approach is that two 
configurations in the update sequence separated by only a few sweeps will still be 
related, or "nearby" in configuration space. So, many sweeps are typically required 
to evolve the sequence between "distant" regions of configuration space. 

As discussed in section 11.41 in lattice QCD with dynamical fermions, the Grass- 
mann degrees of freedom are integrated out analytically, leaving a partition function 
expressed in only in terms of gauge degrees of freedom 



jpf/j g-SG+TrlogMj^-TrlogMpv- (2 4) 

or, equivalently, in terms of gauge and pseudofermion degrees of freedom 

Z = j [VU] [V<^f] e-^G-^^^^^'^^-^^v-*^^^*^^. (2.5) 

Local update processes will no longer be efficient in either case since the calculation of 
A^eflf is no longer local. While solving the sparse linear system MpXf = $f is not as 
computationally demanding as the full diagonalization normally required to compute 
TrlogM, both are essentially non-local computations. So, a successful fermionic 
update process should simultaneously update a large number of degrees of freedom 
for a given calculation of the solution to MfXf = $f without a corresponding large 
change in the action. 

One popular algorithm that accomplishes this is the hybrid Monte Carlo algo- 
rithm. The strategy is to treat the gauge degrees of freedom as generalized coordinates 
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in a Hamiltonian framework of one higher dimension (not to be confused with the 
extra dimension of domain wall fermions). We introduce conjugate momenta H^{x), 
which appear in the extended path integral as a trivial bosonic field 



now serves as the Boltzmann weight for the extended system {U,H). The update 
process for the momenta merely requires the generation of a Gaussian distributed 
field, which is always accepted. Hamilton's equations of motion are then used to 
evolve the system along the classical trajectory. Provided the computational cost of 
the classical evolution is not prohibitive, the result in a simultaneous update of all 
the degrees of freedom between calculations of the (hopefully small) change in the 
Boltzmann weight. 

2.2 The hybrid Monte Carlo algorithm 

The hybrid Monte Carlo (HMC) algorithm combines a non-local update process called 
hybrid molecular dynamics (HMD) with the Metropolis acceptance test described in 
section 12.11 . The full path integral is rewritten in terms of the pseudofermionic 
fields and $py as in (j2.5|l . The algorithm starts with the generation of Gaussian 
distributed momenta H^{x) and the computation of their contribution to the initial 
Hamiltonian. For the fermionic and Pauli-Villars contributions, pseudofermion fields 
$ir and $py are generated for the gauge field configuration using a bosonic heat 
bath algorithm described below. Then, using Hamilton's equations of motion, the 
extended system (f/, H) is evolved in the presence of the forces due to the fixed 
pseudofermion background along the classical trajectory for some finite interval of 




(2.6) 



where the Hamiltonian 




(2.7) 
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the microcanonical evolution parameter, r. Finally, the Hamiltonian is recomputed 
at the end of the trajectory and the change between endpoints ATi is used as input 
to a Metropolis acceptance test, which makes the algorithm exact even though the 
numerical integration only approximately conserves the Hamiltonian. The result is a 
simultaneous update of all dynamical degrees of freedom between each computation 
of the Boltzmann weight. 

For lattice QCD, after formally integrating out the fermionic (and pseudofermionic) 
variables, the only degrees of freedom are the gauge fields. The conjugate momenta 
give the microcanonical evolution of the gauge fields through Hamilton's first equation 
of motion 



The restriction that the gauge links Un{x) remain in the gauge group SU(Ai"c) during 
the evolution constrains the momenta H^{x) to be elements of the Lie algebra of the 
gauge group, hence traceless and Hermitian. Expectation values will be unaffected 
by the shift in normalization of the partition function (|2.6|) due to the inclusion of 
the Gaussian distributed H^{x). 

Hamilton's second equation of motion is derived from the requirement that the 
Hamiltonian is a constant of the motion 



In section 1231 we will explicitly compute d^cfr/dr using ()2.8p and show it can be 
written as 



where -F^(a:) is the called the HMD force. The restriction that the momenta H^{x) 
remain in the Lie algebra of the gauge group during the evolution means that only 
the traceless anti-Hcrmitian part of the force may influence the momenta. Putting 



U^{x) = iHfj_{x)U^{x), (no sum over /i). 



(2.8) 




(2.9) 




(2.10) 
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()2.9|) and ()2.10|) together gives us Hamilton's second equation 

^d^ = (2.11) 

where the subscript TA means to take the traceless anti-Hermitian part. In section 
12. 4[ we will show that each calculation of the fermionic force will require solving the 
sparse linear system MpMpXF = 

The HMD update process starts with the initial gauge configuration Ui and as 
mentioned above, the initial momenta Hi and pseudofermions are generated. This 
completely defines the initial state of the dynamical variables {Ui,Hi) in the back- 
ground environment of $j7 and $py. The gauge and momentum contributions to the 
initial Hamiltonian Tij are computed directly. The pseudofermions are produced by a 
bosonic heat bath, which consists of generating random complex fields Rp and Rpv 
with a Gaussian distribution proportional to exp{—R*R). Defining the pseudofermion 
fields 

<^F = M^Rf, $py = {Mpl^yRpv (2.12) 



-^F (^mI-Mf^ ^ $F 



and 



will give probability distributions proportional to exp 

exp (^—^pyMpyMpv^pv^- Using the Hermitian matrices M^M allow for the use 
of efficient algorithms for solving sparse linear systems. For both Wilson and do- 
main wall fermions, det M^M = det so the algorithm will generate distribu- 
tions for two degenerate flavors of fermions. Note that the computation of $i? is 
local whereas the computation of $py first requires solving the sparse linear system 
(jidpyMpv^ = Rpv and then computing $py = Mpy^'py- The fermionic 

and Pauli-Villars contributions to the initial Hamiltonian are then computed using 
the pseudofermion fields. 

The next step in the update process is to evolve {Ui,Hi) by integrating both of 
Hamilton's equations simultaneously over some finite interval of [rj, Tj+i] to arrive 
deterministically at the final state (f/j+i, i/j+i) ^ {Ui,Hi), all the while holding 
fixed the pseudofermion fields. To satisfy detailed balance, the evolution must be 



40 



reversible: starting from the endpoint of the classical trajectory with the momenta 
reversed, (t/j+i, — i/j+i), and integrating over the interval [Tj_|_i,rj] returns us to the 
starting point, {Ui,—Hi) (tZj+i, — /fj+i). While a continuous trajectory will be 
trivially reversible, on the computer care must be taken to ensure that integration 
using a finite difference approximation is also reversible. 

For the simulations in chapter El the reversible finite difference approximation we 
used was the leap frog method. We choose the trajectory to start at r = and end at 
a fixed r and divide the interval into steps of equal length step size At. To start 
the leap frog, we evolve the momenta one half time step using 

dH (r 0) At 

iH^ix, At/2) = tH^ix, 0) + z ' + 0{At^). (2.13) 

ar z 

We then alternate evolving the gauge field a single time step 

U^{x, t + At)= e'^''("'"+^"/2)Ar^^(^^ ^ 0{At^) (2.14) 

and the momenta a single time step 

r\H (r t) 

iH^ix, T + At/2) = iH^ix, t - At/2) + z ' ^ Ar + 0{At^). (2.15) 

dr 

We finish the trajectory after evolving the gauge field N times with a final half step 
evolution of the momenta. We emphasize that each time we evolve the gauge fields 
we must solve again the linear system MpMpXF = $f and recompute the fermionic 
force needed to evolve the momenta (see section injl . 

To ensure ergodicity, the number of steps must be of order Ar~^, so the differ- 
ence between the initial and final Hamiltonians, ATi, will be of order Ar^. ATi is 
then used in a Metropolis acceptance test at the end of the HMD trajectory to com- 
pensate for finite Ar errors. As the Hamiltonian is an extensive quantity, the step 
size must be decreased as the volume is increased to maintain a given acceptance. 
After the acceptance test, the algorithm continues by generating new momenta and 
pseudofermion fields and evolving along a new classical HMD trajectory. So, as 
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promised, the HMC algorithm can make simultaneous acceptable changes to all the 
gauge degrees of freedom and compute changes to the Boltzmann weight each time 
the fermionic linear system is solved. However, as the lattice volume grows, more 
solutions of the fermionic system are required to produce independent configurations. 

2.3 The preconditioned domain wall action 

Because of the central role of algorithms for solving sparse linear systems in the 
evolution of fermionic configurations, we will reformulate the domain wall fermion 
matrix into a form better suited for more efficient computation. In general, the cost 
of solving linear systems is proportional to the ratio of the largest eigenvalue Amax to 
the smallest eigenvalue Amin of the matrix, called the condition number. A technique, 
called preconditioning, seeks to transform the given matrix (usually at some cost) 
into a new matrix with a smaller condition number. If the savings in solving the 
new system outweigh the cost of transforming the results back into the old system, 
then the preconditioner is successful. For a general discussion on the importance 
of preconditioning when solving linear systems, see Golub and Van Loan [IH]- For 
fermionic evolution algorithms, if the preconditioning transformation preserves the 
determinant, then the cost of the transformation is irrelevant because we can choose 
to work only with the new system. In the case of the HMC algorithm, preconditioning 
will help reduce the computational cost by at least a factor of two. 

The Wilson fermion matrix ()1.3H) contains entries which connect only nearest 
neighbor sites. We introduce a labeling scheme where a lattice site is labeled even 
if all its nearest neighbor sites are labeled odd and vice-versa. For example, if 
^^^j^ x^mod2 = on a given site we label it even (e), otherwise we label it odd 
(o). We then rearrange the components of the fermion vectors to list odd ones first 
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and write the matrix in block form 



loo ^^oe 



M(^) = I I . (2.16) 



Note that we have rescaled the fermion fields by where 

1 



(2.17) 



2(mo + dr) ' 

which makes the on-diagonal blocks the identity. 

DeGrand [13 HHj proposed a preconditioning technique that replaces the Wilson 
fermion matrix M'^'^^ by a preconditioned matrix M^"^^ which preserves the deter- 
minant det M^"^^ = det M'^^^ . He noted that M*^^-' could be decomposed in block 
form 

M(^) = UM^^^L (2.18) 

by choosing U and L to be the upper and lower block-triangular parts of M*^^\ Then 
the block-diagonal M^^) is 

Mm = u-'M^'^^L-' ^ / loo - ^'D^eDpo \ ^^^^g^ 

V Oeo tee J 

We can see that det U = det L = 1 so the transformation preserves the determinant. 

To see that this odd-even technique reduces the condition number, we consider 
an eigenvalue of D, the off-diagonal part of ()2.16|) . where Xd = l/(/t + e). The 
corresponding eigenvalue of M^^^ is Xm = e/{K + e). We can choose e -C k so 
Amin(^) = X^[ ~ e/K. For each Ad, there is a corresponding eigenvalue A|, of 
the block matrix DoeD^o of ()2.19|) . So, the smallest eigenvalue of the preconditioned 
matrix is Amin(M) = 1 — k^/ (/t+e)^ ^ 2e/K and the condition number is approximately 
reduced by half. 

One additional advantage of this preconditioning choice is the decoupling of even 
site fermion degrees of freedom which allows us to integrate them out trivially, adding 
an irrelevant factor to the overall normalization of the partition function. This reduces 
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the storage requirements when implementing the action on a computer, since only half 
of the original components contribute. 

In chapter ^ we relied on our experience with Wilson fermions to guide us in 
constructing and interpreting the physics of domain wall fermions. Here, we can rely 
on that same experience to guide us in rewriting the domain wall fermion matrix in a 
form better suited for solving linear systems. The domain wall fermion matrix dI ^, 
in ()1.42j) will have the same structure as M'^'^^ in ()2.16|) when the fermion fields are 
rescaled by y/R^ where 



Before we continue, it is useful to note that the odd-even block matrices on the d+1 
dimensional lattice will still have a remaining Wilson-like odd-even substructure on 
each d dimensional plane of fixed s coordinate, which follows from ()1.38j) . Since Wil- 
son fermions have been studied by the lattice community for a long time, it is rather 
straight forward to implement an efficient Wilson inverter and verify its implemen- 
tation against the published results of many different groups. So, reusing the same 
code to apply these Wilson Doe matrices to d dimensional sub-vectors in the domain 
wall context will require arranging those components in a separate d dimensional 
odd-even storage order. As a schematic, if the subscripts o, e indicate odd and even 
blocks on the d dimensional lattice and O, E indicate odd and even blocks on the 
d + 1 dimensional lattice then a useful arrangement of vector components is 



1 



(2.20) 





v 



V 




With vector components written in this way, 



of ()1.41|) can be written in 
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block form 



][fidwf) 




(2.22) 



In analogy with ()2.19|) . the preconditioned domain wall fermion matrix, in block 
diagonal form, is 



It is straightforward to devise a scheme to test the implementation of preconditioning. 
After implementing the action of Dqe and Deo on even and odd vectors, one can 
simply verify the following is true: 



for arbitrary \E'. The tests we have performed indicate preconditioning was imple- 
mented within machine precision. For the numerical simulations we present in chap- 
ter |21 the preconditioned domain wall fermion linear solver typically converged in half 
as many iterations as the unpreconditioned one. 

2.4 The hybrid molecular dynamics force term 

In section 12.21 we discussed aspects of the HMC algorithm that did not depend on 
the explicit form of the action. Here we present an explicit derivation of the HMD 
force Ffj^{x), as defined in ()2.10|) . for which we need to compute the derivatives with 
respect to the microcanonical evolution parameter r of all the terms in the action. We 
will consider the Wilson plaquette gauge action, (|1.23p . the plaquette plus rectangle 
gauge action ()1.25|) and the preconditioned domain wall action described in section 
12. 3[ including the Pauli-Villars term. The gauge action contributions to the force are 
well known and are included here only for completeness. 




(2.23) 



(2.24) 
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To determine the gauge contribution to the force, we must compute the derivative 
of the gauge action, which for either action can be written 

dSc P 



dr 



2N, 



E 



tr 



V,{x) + VUx) 



(2.25) 



The Vf^{x) are called staples and are made of sums of SU(A^'c) matrices. We write the 
staple field as the weighted sum of the plaquette and rectangle staple fields: V^{x) = 
coVii^\x) +ciVl;,^\x). For the plaquette field defined in ()1.24p the corresponding 
plaquette staple field is 



[X 



+ Ul{x + fi- v)Ul{x - v)U^{x - v)] 



(2.26) 



For the rectangle field defined in ()1.26|) the corresponding rectangle staple field is 



X 



E 



[U^{x + fi)U^{x + 2fi)Ul{x + /i + i))Ul{x + y)Ul{x) (2.27) 

+U^{x + fi)Ul{x + 2fi- v)Ul{x + fx- v)Ul{x - v)U^{x - z>) 

+Uy{x + fi)Ui,{x + + 0)Ul{x + 2u)Ul{x + y)Ul{x) 

+Ul{x + fi- v)Ul{x + fi- 20)Ul{x - 2v)Uy{x - 2v)U^ (x — u) 

+U^{x + fl)Ul{x + i>)Ul{x -fl + v)Ul{x - fl)U^{x - /t) 

+ Ul{x + fi — z>)f/^(x — z>)f/^(x — fi — i')U^{x — /i — z>)f/^(x — /t)l . 



The expression (j2.25p is valid for either gauge action considered, as the plaquette 
action result is given by setting ci = 0. 

To identify y^), whose traceless anti-Hermitian part will be the pure gauge 
contribution to the HMD force, we substitute Hamilton's first equation of motion 
()2.8p into the derivative of the gauge action (j2.25|) 

dSc P 



dr 2Nr 



J2tH^H,{x) [U,{x)V,{x)-Vi{x)Ul{x)]}. (2.28) 



X,fJ, 
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By comparing this relation with Hamilton's second equation of motion in ()2.10|) and 
(j2.11|) . we choose the Fjf^ to implement in our numerical simulations to be 

= -j^U,{x)V,{x), (2.29) 

which has the same traceless, anti-Hermitian part as the expression in square brackets 
above. We will make similar simplifying choices for the numerical implementation of 
the other HMD forces. 

It will be convenient as we compute the domain wall fermion force term to express 
the result at each step as two terms related by Hermitian conjugation. This will be 
easier if we define 

g{±l2)^^s.x',s' = (1 =F lfi)Sx±fi,x'Ss,s' (2.30) 

and 

Gx,s;x',s' =[!-(! + mf)5o^s'] P+5x,x'Ss+l,s' + ^gi+IJ')x,s;x',s'Uf,{x) (2.31) 

then we can write the preconditioned domain wall fermion matrix which acts only on 
odd sites as 

M = loo- 4DoeDeo (2.32) 

where 

Dx,s;x',s' = Gx,s;x',s' + 75'^s,s"G'^,^„.^, ,,,„757^s"/,s/ (2.33) 

and 7ls,s' = Ss,Ls-i-s' is the reflection operator along the extra direction. Now, it is 
easy to verify the Hermiticity relations = •y^TZDj^Tl and = •y^TZM'y^Tl. 

As with the gauge part, we must compute the derivatives with respect to micro- 
canonical time T of the fermion and Pauli-Villars parts of the domain wall action, 
both expressed in terms of pseudofermions. At first, we work with only the fermionic 
part, Sp, of the domain wall action so we suppress the F and PV subscripts used 
in sections 12.11 and 12.21 until later in this section when we consider the Pauli-Villars 
part. An efficient iterative method is used to solve M'^Mx = in our case based on 
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the conjugate gradient algorithm 
in terms of the solution x as 



The derivative of Sp can then be written 



-X' 



(2.34) 



Using ()2.33p . the derivative of M can be written 
dM 



"d7" 



4 



dr 



+ 



dr 



Deo + I^oi?757^^^757^ + ^^n^^^^KDEo 



dr 



dr 



If we define another quantity 



AGeo . '^GoE 

^OE 1 J-^EO 



dr 



dr 



(2.35) 



(2.36) 



then the derivative ()2.35|) can be rewritten 

dM 



dr 



(2.37) 



after using the Hermiticity relations and noting {^^TVj = U. From here, it is not 
difficult to show that 



dr 



757^/^757^M + j^lZM-f^nK + h.c. 



X, 



(2.38) 



accomplishing our goal of expressing the derivative as two terms related by Hermitian 
conjugation. 



The next step is to arrange the derivative calculations using ()2.8|) so that iH^ 



x] 



appears on the left. We start again from ()2.31|) 

-^Gx,s;x',s' = 'i'H^j,{x)g{+lj)x,s;x',s'Ufi 



X]. 



(2.39) 



While G depends explicitly on the bare mass m^, the derivative does not as there 
are no gauge fields in the extra direction. This result is then substituted into ()2.36p . 
recalling that K is defined only when both sets of indices are odd 



x,s;x's' 



-t^l X] [Dx,s;x",s"iH^{x")g{+^)x",s";x',s'U,,{x") 



+iH^{x)g{+fl)x,s-x",s"U^{x)Dx",s";x',s 



(2.40) 
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Finally, we substitute into the derivative of the fermionic action, and use the invari- 
ance of the trace under cyclic permutation to move the conjugate momenta to the 
left 

dSp 2\ ^ / ■ TT / \ 

x,n 

X {U^{x) tr,,spin [gi+p) (xX^MW + Dxx^M^) 

+ g\-ii) (Mxx^D^ + D^Mxx^)] (2-41) 

-h.c.}). 

As with the pure gauge contribution, we compare this relation with Hamilton's second 
equation in (|2.1(Jj) and (|2.1H) and choose Flf\x) to have the same traceless, anti- 
Hermitian part as the expression in the curly braces but in a form convenient for 
numerical simulation. 

We construct some vectors which are defined on half of the sites of the lattice 

PE = DeoXo (2.42) 
i}o= -nlMooXo (2.43) 
aE= D^oi^o = -kID^eoMooXo. (2.44) 

Now we can write the fermionic contribution to the force 

Flf\x) = -2U,{x) tr,,3pin b(+/i) (X^^ + P^^) + gK-^^) (V + crx^)] ■ (2.45) 

We can make one further simplification by defining two vectors which extend over the 
whole lattice 

(2.46) 

\i'o J 

Finally, in terms of the vectors v and w 

Ff )(x) = -2U^{x) J2 trspin [(1 - 7,)v,+^,swls + (1 + 7m)^x+a,.^^L] • (2.47) 

s 

With the absence of gauge fields in the s-direction, it is not surprising that this 
fermionic force is identical in form to the force due to flavors of Wilson fermions. 
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This implies that all of the special properties of domain wall fermions, at this level, 
are encoded at the start when finding the solution xo of the fermionic system. 

The derivative of the Pauli-Villars part of the domain wall action is 



dr 



PV 



$py. (2.48) 



Note the crucial difference in sign relative to the derivative of the fermion action 
fl2.34|l . The construction of the Pauli-Villars contribution to the HMD force mirrors 
the construction for the fermionic part provided we set the solution xpv to be the 
pseudofermion field generated by the heat bath, xpv = ^pv, set rnf=l and keep 
track of the overall change in sign. Using xpv in place of x we get 

Ff^)(x) = -Ff)(x)| ^ . (2.49) 
Now that all three force contributions are known, 1)2.111) can be written 

^^^.(^) = [^f n^) +^r(^)]TA- (2-50) 

We again emphasize that only the traceless, anti-Hermitian part of F^[x) contributes 
to the molecular dynamics evolution. 

With any large coding project, it is important to devise testing strategies to verify 
that the code is consistent with the analytic description. For the dynamical domain 
wall action, this process is quite simple and straight forward, provided one starts with 
an accurate implementation of the Wilson fermion matrix acting on a vector, some 
kind of linear solver, and the Wilson HMD force term. Before we implemented the 
domain wall action in the physics software environment of the QCDSP machines at 
Columbia, we ran several dynamical Wilson fermion simulations which were in good 
statistical agreement with the published results of other groups. 

Since our linear solver, based on the conjugate gradient (CG) algorithm, solves 
the problem M'^Mx = testing it simply requires choosing many different $ sources. 
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solving each one for the corresponding x solutions and then applying M to x ^^id 
comparing the result to the original. This defines a quantity called the true residual 

$ - WMx 

(2.51) 



|<l>| 

where we use the vector 2-norm: |$| = V However, running the CG algorithm 
for a sufficient number of iterations for r to be smaller than some predetermined 
stopping condition only verifies that some non-singular linear system was solved. 

To check that the linear system solved was actually involved the domain wall 
fermion matrix Mt(dwf)]Vf (<i™f) acting on a vector, we compared the results of the 
solution run on a free lattice configuration, Un{x) ^ 1, to the analytic solution com- 
puted by Vranas 14 . We then extended our check to non-trivial gauge configurations 
by applying a random gauge transformation, A(x) G SU(Ai"c), to the free lattice con- 
figuration, Ufi{x) A'^{x)A{x + fi), and verified that the true residual was invariant 
under gauge transformation 

(2.52) 



|A$o| 

where $0; Xo ^^e the free field source and solution verified by analytic calculation 
and Ma indicates the matrix is applied using the gauge transformed links. If the 
domain wall fermion matrix is properly implemented, then the true residual is a 
gauge-invariant quantity. When starting with a properly implemented Wilson fermion 
matrix, the gauge invariance test should not fail if the free field test passed, since the 
gauge fields only enter the domain wall matrix through the Wilson portion of the 
code. These tests passed within finite precision limits of the QCDSP machine. 

The only difference between the force terms of the fermionic ()2.47|) and Pauli- 
Villars (j2.49|) parts of the domain wall action and the force term of Lg independent 
fiavors of Wilson fermions is in the computation of the input solution x- So, start- 
ing from an accurate implementation of dynamical Wilson fermions that supports 
multiple independent fiavors makes the implementation of dynamical domain wall 
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fermions straightforward. In essence, once x has been computed by the hnear solver, 
s is treated as a Wilson flavor index. As a final check, we also verified that the domain 
wall force terms transform properly under random gauge transformations. 

After these series of analytic tests of the numerical implementation, it would have 
been nice to run an actual simulation and do a statistical comparison with published 
results. Since the simulations described here are the first of domain wall fermions 
in dynamical QCD, there are no previous published results. However, there were 
still two significant statistical tests which were run. First, we simulated dynamical 
domain wall fermions at fnf=l on 8'^x32 volumes at [3=5.7 and verified the results 
were consistent with quenched QCD. This check is not as trivial as it may seem, as 
the two parts of the domain wall force do not exactly cancel trajectory by trajectory 
since different pseudofermion fields are used for the fermionic and Pauli-Villars parts. 
Instead, the cancellation of forces happens stochastically, yielding quenched results 
in the limit of a large number of trajectories. Second, Vranas performed dynamical 
domain wall simulations of Nf=2 QCD on 2^ volumes 50^ and found the results were 
in good statistical agreement with dynamical simulations done directly in the overlap 
formalism. 

Perhaps the best test of the whole implementation of the CG solver and domain 
wall force working together within the context of the HMC algorithm comes from an 
application of Liouville's theorem by Creutz [3] within the classical framework of 
molecular dynamics. If we integrate out all of the fermion (or pseudofermion) fields, 
we can write the molecular dynamics partition function as a path integral over gauge 
fields and conjugate momenta 

ZuMD = J [VW] [VH'] e-^' = J [VU] [VH] e-^e«-^'. (2.53) 

Since the HMD evolution is purely classical, the area preserving property of Liou- 
ville's theorem guarantees that the integration measures are equal: [DW] [DH'] = 
[VU] [DH]. Dividing both sides by the partition function, Zhmd, we get the expec- 
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tat ion values 



(exp(-AH)) = 1 



{AH = n'-n) 



(2.54) 



and, by Jensen's inequality, 



{AH) > 0. 



(2.55) 



Because ATi, is computed as part of the Metropolis acceptance test of the HMC 
algorithm, it should be straightforward to calculate these observables in any actual 
simulation. 

In our experience, satisfying these identities seems to be the most stringent check 
of the implementation of the HMC algorithm. If the force term is inconsistent with 
the Hamiltonian as implemented then Liouville's theorem will not hold because the 
evolution is no longer classical and area preserving with respect to the given Hamil- 
tonian. However, the code necessary to compute the change in the Hamiltonian need 
not share any components with the force term except for the linear solver. Once the 
CG solver has been sufficiently tested, checking the Hamiltonian part of the code is 
simple, so these identities provide a demanding statistical test of the overall evolution, 
and in particular the HMD forces. 

In chapter El the description of our simulations will include the observable ()2.54|) . 
All of the runs, taken as a whole, agree with the prediction ()2.54|) within statistical 
error, from which we conclude that the domain wall HMC algorithm was correctly 
implemented. 
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Chapter 3 

Simulations of Nf=2 QCD 

3.1 The deconfinement phase transition 

We are motivated to study the thermodynamics of Nf=2 QCD with degenerate hght 
quarks for two reasons. First, it is an approximation to Nf=6 QCD with the six quarks 
at their physical masses. Recall from section lL6l the discussion of the quenched limit 
of QCD, where dynamical quark effects are altogether eliminated by taking the bare 
quark masses to infinity. Similarly, the dynamics of the heavier quarks {top, bottom, 
■ ■ ■) will be marginal in the Euclidean path integral and hence will have little effect 
on the character of the transition region. So, Nf=2 QCD is the most natural starting 
point for studying the role of dynamical quark effects in QCD thermodynamics. 

At first glance, this theory, corresponding to a world where the up and down quark 
masses are degenerate and the rest are much heavier, is a reasonable approximation 
of our world. However, on closer inspection, the up and down quark masses are not 
exactly degenerate and the strange quark mass is only an order of magnitude heavier. 
So it will be important to perform the typically more computationally demanding 
Nf=3 QCD simulations with non-degenerate quarks to understand the full physical 
theory in detail. But, one should first study completely (if possible) the simpler Nf=2 
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theory before looking for (perhaps) subtle differences by turning on the dynamics of 
the heavier quark flavors. 

A second motivation for this study is simply that Nf=2 QCD is an interesting 
theory in its own right. With massless quarks, the theory minimally has 811^(2) ® 
SU/j(2) chiral symmetry and the chiral condensate, (qq), is an exact order parameter. 
This chiral symmetry group is expected to be in the 0(4) universality class and other 
models in that class are known to have second order phase transitions, so there is a 
possibility that the chiral transition of Nf=2 massless QCD may be second order |28j. 
However, if the anomalously broken U^(l) symmetry is effectively restored along with 
chiral symmetry at the transition, then the larger symmetry is expected to be in the 
0(2)00(4) universality class and as other models in that class are known to have 
first order phase transitions there is the possibility that the chiral transition could be 
first order [2HlEniEIl- 

Before embarking on a large numerical project with a new fermion formulation, we 
note what is already known about QCD thermodynamics using other formulations, as 
all fermion formulations should give the same answer in the continuum limit. On 24^x4 
volumes in quenched QCD, there is a first order deconfinement phase transition at 
/?c=5.6925 • Along with other quenched simulations, this corresponds in physical 
units to a critical temperature Tc~ 270 MeV [SB]- 

For Nf=2 QCD with staggered fermions on 16^x4 volumes, the finite temperature 
phase transition occurs at /5c=5.265 for m=0.01 and /?c=5-291 for m=0.025 O EH] 
which suggests that /3c ~ 5.25 in the chiral limit, approximately 150 MeV. At these 
quark masses, there is no evidence for a first order phase transition. Other Nf=2 
dynamical simulations with different actions have critical temperatures in physical 
units in the range ~ (170-190) MeV 

At finite lattice spacing the chiral symmetry of staggered fermions is only a one 
dimensional subgroup of the continuum symmetry and the vector flavor symmetry is 
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explicitly broken. Consequently, QCD with staggered fermions at accessible lattice 
spacings is a theory with one physically light pion and two unphysically heavy pions 
in the range 500-600 MeV. This may obscure the true nature of the transition until 
the lattice spacing is reduced or the action is improved. 

In section we will locate the transition region of Nf=2 QCD with moderately 
light degenerate quarks in small 8^ x 4 volumes. Using the chiral condensate we 
will show strong evidence for a chiral phase transition in the massless quark limit. 
However, the small volumes used will make it difficult to establish detailed properties 
of the transition region. Thus, in section IT^ we will focus in on the transition region 
in larger 16^ x 4 volumes with somewhat lighter quarks in the hope of seeing an 
indication of the order of the transition. 

Implicit throughout the preceding discussion is the assumption that the domain 
wall fermion formulation provides us with the nearly exact chiral symmetry needed to 
observe the critical phenomena related to spontaneous chiral symmetry breaking. In 
section Em we will perform the analysis of the GMOR relation outlined in section fT31 
as a quantitative measure of the effectiveness of the domain wall method for producing 
nearly massless fermions on the lattice. 

3.2 Locating the transition in 8^x4 volumes 

In the Euclidean path integral formulation, the temporal extent of the lattice Lf is re- 
lated to the inverse temperature through the lattice spacing: Lt = (aT)^^. If we 
hold Lt fixed, then the determination of the temperature Tc of the chiral phase transi- 
tion requires, in general, searching in the space of the four parameters (/?, mo, Ls,mf) 
that determine the lattice spacing. The problem of searching this space is made sim- 
pler as we are only interested in the location of the transition in the chiral limit, as 
ruf ^ and Lg —>■ oo, over an appropriate range of mo. In the massless limit, the 
location is specified by the critical coupling, Pd^o)- 
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We first studied the behavior of both the Wilson hne and chiral condensate in 
8^ X 4 volumes gUl iH iSl 1121 EEI EO] • Choosing the spatial extent of the lattice L to 
be much larger than the temporal extent, L^Lt, is essential to eliminating the finite 
volume effects that obscure the character of the transition. However, past experience 
and this work will demonstrate that searching in small volumes with moderately heavy 
quarks gives a surprisingly good estimate of f3c- Further simulations at larger spatial 
volumes and smaller masses can refine the estimate of /?c and reveal the character of 
the transition. 

At very small lattice spacings, when the gauge fields are near the identity, the 
condition for bound domain wall states is essentially the free condition ()1.49|) . For 
a large range of mo, one light flavor will be localized on the domain wall defect. 
Of course, to use the HMC algorithm described in section 12. 2| we must square the 
determinant which gives us two degenerate flavors. However, mo also determines 
the maximum four-momentum a bound surface state may have, see ()1.49|) for the 
condition in the free case, and is the mass scale of the heavy Wilson doubler states 
that are not localized on the defect. So, at coarser lattice spacing, these scales become 
relevant and the bare couplings and masses will have mo-dependent contributions 
which should cause shifts in (3c- By locating Pc for several mo, we hope to establish 
that this parameter will not have to be fine-tuned on coarse lattices for the domain 
wall method to work, but some care must be taken in choosing a value to ensure 
the simulation represents the correct number of light flavors with a sufficiently large 
phase space for the light degrees of freedom. 

Consideration of the existing numerical studies of QCD thermodynamics with 
quenched and dynamical Nf=2 staggered fermions provides us with a rough esti- 
mate of the location of (3c- For quenched Lf = 4 thermodynamics, corresponding to 
rrif = 1 in our parameter space, /3c ~ 5.7 and for dynamical staggered fermions in 
the chiral limit, similar to mj = and Lg — > oo, (3c ~ 5.25. Since we are interested in 
the thermodynamics of domain wall fermions near the chiral limit, if we choose rrif 
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sufficiently small and Lg sufficiently large then for some reasonable range of mo we 
should expect to find /3c in the vicinity of the chiral staggered value. 

To perform a fully dynamical search of the parameter space requires starting a 
new simulation each time any one of the basic run parameters {L, Lt, l3,mo, Ls,mf) 
is changed. Although the process of analyzing any given simulation is interesting in 
its own right, the sheer number of simulations included in this work prevent us from 
providing that level of detail. So, we encourage the reader to refer to the extensive 
tables and figures provided for details of any particular evolution. 

Our search strategy for 8'^x4 volumes was to fix ruf =0.1, Ls=12 and for several 
values of itlq from 1.15 to 2.4, we scan in j3 from weak to strong coupling over the range 
j3 = 5.95 - 4.65 until crossover behavior was observed in (qq) and (11^31). The run 
parameters for the simulations are listed in tables^-IHl There were no difficulties with 
the CG algorithm failing to converge, characteristic of the "exceptional configuration" 
problem in dynamical Wilson fermions. When the even-odd preconditioned domain 
wall Dirac operator was used, as described in section 12.31 100-130 CG iterations 
were typically needed per HMD step in the confined phase and 40-70 were needed 
in the deconfined phase. Simulations took three to five days to complete 800 HMC 
trajectories when run on a six Gflops portion of the QCDSP computer at Columbia. 
In some cases, when preconditioning was not used, two and a half times the number 
of CG iterations were needed per HMD step. 

As the deconfined phase exists at weaker coupling than the confined phase, the 
gauge configurations are more ordered, i. e. closer to the identity configuration. So, 
the initial configuration for each simulation was deliberately chosen to be as far as 
possible from the anticipated state after thermalization. For simulations where a 
deconfined thermal state was expected, the initial lattice was disordered, while an 
initial ordered lattice was used where a confined thermal state was expected. This way, 
the thermalization time scale could be estimated by observing the system tunneling 
from the wrong phase through the transition into the correct phase. 
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The evolutions of (gg) and (iH^sl) are shown in figures ^-^1 The vertical dashed 
line indicates that data points to the left are considered part of thermalization and 
were not used to compute averages and autocorrelated errors. The figures clearly show 
that the short (microcanonical) time fiuctuations are much larger in the confined 
phase than in the deconfined phase. Note that the data show a tendency for the 
thermalization time and the scale of large time fiuctuations to increase towards the 
crossover region. This is an encouraging sign that correlation lengths are increasing 
towards the center of the crossover region, suggesting the long range order typical of 
phase transitions. 

The averages and autocorrelated errors of (gg) and (IW3I) are plotted versus [3 
for mo = 1.15 - 1.8 in figure IT7I and mo = 1.9 - 2.4 in figure ITHl For mo > 1.65 
there is a clear indication of simultaneous rapid change in (gg) and (IM/3I) versus (3, 
indicative of a potential phase transition smoothed out by the moderately light quark 
mass and small volume. The reader should be careful to note the change in scale for 
(gg) between the two figures. For mo < 1.4, there is still a crossover signal in (IVF3I), 
suspiciously close to the transition region for quenched QCD, but any signal in (gg) 
is nearly gone. 

Once we have identified the crossover region, we can estimate the location of j3c 
through a fit of the approximate order parameters to the ansatz 

f{(3) = Co {ci + tanh [c2 {(3 - (3,)]] . (3.1) 

The results of the fits are in tables IHl and QHl plots of the fit curves versus the data 
are also shown in figures [T7| and UHl and a plot of [3c versus mo is in figure UHl 

Although the purpose of the ansatz ()3.1|) is to estimate the infiection point /3c, 
the other parameters of the fit can be used to quantitatively describe what our eyes 
have already told us by looking at the figures. After extracting the overall scale factor 
Co we can compare ci and C2 between the various fits. When sign(co)ci = 1 then the 
fit describes an approximate order parameter that goes exactly to zero deep in the 
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the symmetric phase. C0C2 is the slope of the fit curve at the infiection point, so C2 
indicates the relative rate of crossover between the two phases. Examining table El 
we see from c\ that the Wilson line is a good approximate order parameter and from 
C2 that there is a tendency for the transition to broaden as mg increases. 

In table ^1 we see more interesting behavior. First, we note in figure El that 
the /3-dependence of (gg) for mo < 1.4 is not well modeled by The tendency 

of I Co I to increase with increasing mo is expected since this is just the change in the 
overlap between our effective quark fields ()1.51|) and the actual light fermion modes 
suggested by ()1.50|1 . The vanishing of this overlap as mo < 1.4 is our first indication 
that any light fermion modes are, on average, no longer localized to the domain wall 
defect. Next, we note that |cl| ^ 1 for mo < 1.4 which indicates that (gg) is no 
longer approaching zero in the symmetric phase and again suggests the phase space 
for light fermions on the defect is nearly gone. We conclude that domain wall fermions 
are unreliable for Nf=2 and Lt=A thermodynamics with mo < 1.4. Because the fit 
ansatz is only intended to give a qualitative description of the data, we expect that 
X^/^dof ^ 10 indicates a reasonable agreement between the data and the ansatz. 

Finally, for C2 we see the same tendency for the transition to broaden towards 
larger mo. We expect that for some mo > 2.0 the simulations should behave like Nf=8 
QCD but our data do not indicate any sharp crossover from Nf=2 QCD to Nf=8 
QCD for mo < 2.4. One interpretation is that the simulations are Nf=2 with the shift 
in jSc due to the change in the amount of phase space available to the light fermions 
as mo varies, in analogy to the free field condition ()1.49|) . A related interpretation is 
that as the phase space for the two light modes decreases while the modes themselves 
become heavy for mo > 2.0 the phase space for the eight doubler modes is increasing 
while the doubler modes become light, so that mo smoothly interpolates between two 
fiavor and eight fiavor thermodynamics. 

Although the rapid crossovers seen in 8^x4 volumes with Ls=12 and m/=0.1 are 
indicative of a phase transition, the dynamical quarks are likely too heavy to expect 
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good chiral behavior. Furthermore, it is difficult at this point to quantify the residual 
chiral symmetry breaking effects due to finite Lg. We can strengthen our evidence for 
the existence of the phase transition by performing a dynamical extrapolation to the 
massless quark limit: oo, mj — > 0. Thus, if we study the order parameter (gg) 

outside the transition region and demonstrate that it is zero in the symmetric phase 
and non-zero in the broken phase then we are reassured that the rapid cross over is 
strong evidence for a phase transition in the chiral limit. 

We choose mo=1.9 to study the separate phases for several reasons. First, mo=1.9 
is sufficiently far from the region mo < 1.4 where the light surface states of domain 
wall fermions vanishes. Second, choosing mQ\2 guarantees that at low momentum only 
a single light surface mode is bound to the domain wall so ensuring the low momentum 
phase space of an Nf=2 theory after squaring the determinant. Finally, although not 
strictly necessary, the value of f3c is similar to the value for Nf=2 staggered fermions 

To study the confining (broken) phase in an 8^ x 4 volume with mo=1.9, we 
choose (3=5.2 as a good compromise between the desire to work at as weak a coupling 
as possible to suppress the residual chiral chiral symmetry breaking effects of finite 
Lg and the need to be far enough from the crossover region that the results will 
be unaffected by the presence of a nearby phase transition. To perform the chiral 
extrapolation, we varied the bare quark mass ruf over the range 0.02 to 0.18 and the 
extent of the extra dimension Lg over the range 4 to 40. The run parameters and 
results are summarized in table El 

Focusing on the chiral condensate, we expect from perturbation theory that (qq) 
will be linear in the dynamical quark mass, rrif. In figure EDI we show constrained 
ffis to both a linear ansatz and a quadratic ansatz for Lg = 8, 16 and fnf = 0.02 - 
0.18. By examining the x^/^dof for each constrained ffi in table El we can see that 
the data tend to be more linear in mj as Lg increases. 
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For rrif = 0.02,0.1 there are sufficient data to extrapolate (qq) to — > oo limit 
using an exponential ansatz. The results of these ffis are shown in table El and figure 
1^ We can see that the fit of the data to the ansatz is reasonable. One might worry 
that the x^/A^dof are somewhat large but it is likely that the true autocorrelated errors 
are underestimated due to the rather short length of evolutions. Using the values of 
(qq) as mj ^0 from the fits in table we can also extrapolate to oo and we 

find 

lim lim {qq) (p = 5.2) = 0.0061(1). (3.2) 

Ls— ►oo my— >0 

Checking to see if the order of limits changes the extrapolated value we take the 
Lg ^ oo limits for ruf = 0.02,0.1 and do an unconstrained linear extrapolation to 
iTLf and find 

lim lim {qq) {p = 5.2) = 0.0061(1). (3.3) 

Thus, we conclude that there is a broken phase for j3 < 5.2 indicated by the {qq) 
order parameter in the chiral limit. 

To study the deconfined (symmetric) phase in 8^x4 volumes with mo=1.9, we 
choose to simulate with (3=5.45. As before, we desire to run at as weak a coupling as 
is reasonable, but it is important not to go too deep into the symmetric phase as our 
physical volume will become too small and too hot to yield useful results. As before 
we varied rrif over 0.02 - 0.18 and Ls over 4 - 32. The run parameters and results 
are summarized in table El In figure (221 we show constrained fits for a linear and 
quadratic ansatz for Ls=8, mj = 0.02 - 0.18 and a fit to a linear ansatz for Ls=lQ 
over the same range of rrif. We see immediately that the data appear much more 
linear than in the broken phase. This may occur because the rate of exponential 
suppression of chiral mixing of the domain wall states increases at weaker coupling. 
Another possibility is that the rate of exponential suppression may be larger in the 
symmetric phase due to the smoother gauge fields. This remains a question open for 
further study. Results of the fits are in table [13 where we have also included fits of 
the Ls=16 data constrained be zero at mf=0. 
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As in the (3 = 5.2 study, we also fit exponential extrapolations to Ls —>■ oo, see 
table El and following the same double limit procedure above, if we do not constrain 
the — > oo, m/ — > limit of {qq) to be zero, we get the extrapolated values 



lim 

-Ls— >oo 

lim 

mf — »0 



lim (qq) {(3 = 5.45) 



lim (qq) {(3 = 5.45) 



0.00010(5). 
0.00011(5). 



(3.4) 
(3.5) 



Furthermore, if we constrain the Lg oo, limit of (qq) to be zero, the quality 

of the fits does not change significantly. Thus, we conclude that our observations are 
explained by a deconfining chiral phase transition with (qq) as the order parameter 
for that transition in the chiral limit. 



3.3 The transition region in 16^x4 volumes 

In section 13. 2^ we provided strong evidence of a chiral phase transition for Nf=2 
domain wall fermions in the massless limit: mj— >0, L^^oo. After having bracketed 
the transition region in smaller spatial volumes, we would like to characterize the 
order of the phase transition. Working in large spatial volumes with light quarks is 
essential to observing the long range correlations typical of true critical behavior. We 
chose to increase the volume from 8^x4 to 16^x4. To make the quarks lighter, we 
chose a larger Ls=2A and a smaller m/=0.02. 

Recalling from QCD with staggered fermions that (3c decreased slightly as the 
quark mass is decreased, we first ran simulations in 8'^ x 4 volumes to map out the 
shift in f3c- Further, having simulations in the transition region in two different spatial 
volumes will allow us to study systematic errors due to finite spatial volumes. The 
simulations were run at /? = 5.2, 5.275, 5.325 and 5.45 with mo=1.9. Table fTTI contains 
the other run parameters for these simulations. As before, there were no difficulties 
with convergence of the CG algorithm. 300-450 CG iterations were typically needed 
per HMD step in the broken phase and transition region and 100-120 CG iterations 
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were needed in the symmetric phase. Simulations in the broken phase took about 
thirty days to complete 800 HMC trajectories when run on a six Gflops partition of 
the QCDSP computer. 

The evolutions of (gg) and (IM/3I) are shown in figures |21] and |2H1 The vertical 
dashed line indicates that data to the left are considered part of thermalization and 
were not used to compute averages and autocorrelated error estimates. Even more 
dramatically than in section 13.21 the short time fluctuations of (gg) are much larger 
in the broken phase than the symmetric phase and the longer time correlations are 
clearly visible at /3=5.325. In figure QUI we show (gg) and (IVF3I) vs. (3 and we see that 
/3=5.325 is clearly in the middle of the transition region. 

For 16^x4 volumes, we chose (3 = 5.25 - 5.35, spaced at 0.025 intervals. Table 
irrn contains the other run parameters. Both ordered and disordered starts were run 
for each set of run parameters until observables measured on last three hundred 
trajectories of each simulation were in good statistical agreement. There were no 
difficulties with CG convergence. 260-300 CG iterations were typically needed per 
HMD step, with more iterations required in the middle of the transition region {(3 ~ 
5.325) than on the edges. To complete 800 HMC trajectories on a 25 Gflops partition 
of the QCDSP computer would take between 48 and 72 days. 

The evolutions of (gg) and (IM/3I) are shown in figures HTl and OHl Since two evo- 
lutions appear in each panel of the figure, separate vertical lines are used to indicate 
the start of the last 300 trajectories for each evolution in the panel. The evolutions 
show no convincing evidence for a two-state signal indicative of a first order phase 
transition. 

The summary plots which display (gg) and (IVF3I) versus (3 are shown in figures 
ESI and 1^ The most obvious feature of the transition for domain wall fermions is that 
it does not appear particularly steep. For comparison, we have also plotted Nf=2 
staggered fermion data ^55j. Based on the different symmetries of the two fermion 
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formulations, the staggered data show the transition when one pion is hght and two 
are heavy whereas the DWF data show the transition for three degenerate pions. 

Comparing the 8^ and 16^ volumes, we also see that the larger volume did reduce 
the short distance fluctuation in the signal. However, the point in the middle of the 
transition region, (3 = 5.325, does not show signs of the critical fluctuations indicative 
of a second order transition nor does it show the two state signal of a first order phase 
transition. Based on the broadness of the transition region and the absence of critical 
fluctuations, we believe that the domain wall states are still too heavy to reveal any 
critical behavior present in the massless theory. 

To bolster our claim, we performed an 8^x32 zero temperature scale setting calcu- 
lation, with /5=5.325, mo=1.9, L,=24, m/=0.02 |SH]. We found that mp=1.18(3) and 
m7r=0. 654(3) in lattice units. Using rrip to set the lattice scale, /3c=5.325 corresponds 
to Tc = 163(4) MeV and = 427(11) MeV in physical units. By comparison, for 
staggered fermions near for Lt = 4, one pion has a mass of m,r ~ 230 MeV and the 
other two are nij,.^ ~ 600 MeV j2Zl ESI- This leads to several possible explanations 
for the difference in the widths of the crossover region for staggered and domain wall 
fermions. One possibility is that the width may depend on the difference in the pion 
spectrum: one light and two heavy pions for staggered vs. three moderately heavy 
pions for domain wall fermions. 

We scanned through the domain wall crossover region by changing f3 holding other 
parameters, particularly rrif and Lg, fixed. We expect that additive renormalizations 
of the bare quark mass, mres in the discussion of section II. 5| will be smaller at 
weaker coupling. The mass of the domain wall pions suggest mres ^ mj so it is 
possible that the effective bare quark mass, = mj + mj-es may also be rapidly 
changing through the transition region, leading to qualitatively different behavior 
than staggered fermions, where the effective bare mass did not change through the 
transition region. This will be discussed further in section 13.41 
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Besides the relatively heavy domain wall quark masses to account for the broad- 
ening of the crossover region relative to staggered fermions, a final possibility is that 
the staggered lattice spacing evolves rapidly as a function of bare coupling in the tran- 
sition region, implying the crossover behavior for domain wall and staggered would 
be similar if shown in physical units. The dependence of the staggered lattice spacing 
on the bare parameters in the transition region was extensively studied by Blum, 
et al. [Snj. Summary plots which display the staggered chiral condensate (xx) ^-^id 
Wilson line (IVF3I) vs. temperature T in MeV are shown in figures l35l and l36l 

As we ran only a single zero temperature domain wall scale setting calculation 
in the transition region, there are insufficient data to make a non-perturbative deter- 
mination of the functional dependence of the lattice spacing on the bare parameters, 
as was done for staggered fermions. Close to the continuum limit, it would only be 
necessary to determine the lattice spacing for one set of bare parameters as the lattice 
spacings at other bare parameters are given by the asymptotic scaling relation for two 
flavors of massless fermions 

oAlat = (87r2/3/29)'^'/'^' exp {-A'k^[3/29) (3.6) 

where Alat is the lattice scale parameter which must be determined by a scale setting 
calculation. The relation ()3.6p is the result of a two loop perturbative calculation and 
only contains terms which are independent of renormalization scheme. However, the 
value of the parameter A will depend upon the renormalization scheme 

Solving ()3.6|) for Alat in our 8^ x 32 scale setting calculation at /3=5.325 men- 
tioned above, we find Alat=1-388(2) MeV^. Holding Alat fixed, we can estimate the 
temperature at other couplings holding the other bare parameters fixed. Although 
the effective bare quark mass will change as we change the coupling, as noted above, 

we emphasize that irip is typically less sensitive to small changes in the quark mass 
^The difference between this Alat and the more famihar ~ 200 MeV is merely due to large 
scheme-dependent factors. 
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than rriTT so we expect that this may not be a large source of systematic error. Using 
this method, we also plot the domain wall chiral condensate (gg) and Wilson line 
(IW3I) vs. temperature T in MeV in figures 1^31 and IH^ Plotting both the staggered 
and domain wall data vs. temperature has not resolved the difference between the 
widths of the crossover regions for the two formulations, assuming that the violations 
of the scaling relation (j3.6p are not large due to the relatively strong couplings. We 
will examine this issue below where we have two scale setting calculations at different 
couplings to estimate the size of the scaling violations. 

After realizing that the pion mass was too heavy to reveal critical phenomena, 
we could, in principle, merely increase Lg until we achieve the desired pion mass. 
Unfortunately, the 16^ x 4, Lg = 24 calculation is essentially the largest scale study 
we were able to conduct in a reasonable amount of time. So, we decided to pursue 
improvements to our lattice action which we hoped would reduce the pion mass 
significantly for roughly the same computational effort. 

Encouraged by published quenched results which indicated that replacing the 
standard Wilson gauge action ()1.23|) with the improved rectangle action ()1.25|) can 
reduce the quenched pion mass [SIIEH], we located the Nf=2, Lt=A transition region 
for this action in 8^ x 4 volumes with ci=-0.331, mo=1.9, L,=24, m/=0.02 [I^IHKj . 
The choice of Ci = —0.331 stems from early work by Iwasaki jnUEl]- Although the 
rationale of his approach probably has little relevance to this work, choosing this 
value eliminates the need for costly study of the effects of varying ci. Simulating with 
the additional term in the gauge action did not significantly increase the computa- 
tional cost relative to the standard action because the cost is dominated by the CG 
algorithm. Table El contains the other run parameters for these simulations. 

The evolutions of (gg) and (IW3I) are shown in figures QUI and IHUl As before, data 
to the left of the vertical dashed line in each panel were considered thermalization 
and not used to compute averages and autocorrelated errors. Ordered and disordered 
starts were run at each (3 until the averages of observables agreed within errors for at 
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least the last 100 trajectories. As with domain wall fermions with the Wilson gauge 
action, the short distance fluctuations of (gg) are diminished in the symmetric phase. 

Following the 8^ x 4 study to locate the transition region, we performed a detailed 
study of the transition region in 16^ x 4 volumes. Tablel^ contains the run parameters 
for these simulations. The evolutions of (gg) and (IW3I) are shown in figures 1^ and 
1221 The summary plot which display (gg) versus (3 is shown in figure EZI for both 
8^ X 4 and 16^ x 4 volumes. Contrary to the observed decrease in the pion mass for 
quenched domain wall fermions with the Iwasaki gauge action at zero temperature, the 
broadening the transition region at finite temperature again indicates a moderately 
heavy pion for dynamical domain wall fermions with Iwasaki gauge action (ci=-0.331) 
as was the case with the Wilson gauge action (ci=0). 

We performed 8^ x 32 zero temperature scale setting calculations in the transition 
region, with /5=1.9 and 2.0, Ci=-0.331, mo=1.9, L,=24 and m/=0.02 [58;. For (3=1.9, 
which appears to be in the central part of the transition region, the dimensionless 
lattice meson masses are mp=1.16(2) and m^=0. 604(3) and, using mp to determine 
the lattice spacing give the physical temperature T= 166(3) MeV and physical pion 
mass m^(f^^^^= 400(7) MeV. This confirms our suspicion that the improved gauge 
action ()1.25|1 does not significantly improve the chiral behavior of dynamical domain 
wall fermions at strong coupling. 

We used the same asymptotic scaling technique as before to extrapolate the lattice 
spacing as a function of the bare coupling (3 and find Alat=98.0(1) MeV. We can 
now directly compare the transition regions of the two domain wall QCD actions. In 
figures EH and EHl we show the chiral condensate (gg) and the Wilson line ( | PI/3 1 ) for 
the standard Wilson gauge action (ci=0) and the Iwasaki gauge action (ci=-0.331). 
While the relative offset between the two curves may be due to uncertainties in the 
overall temperature scale of the scaling functions for each data set, superficially it 
appears that the rectangle improvement in the gauge action did not reduce the width 
of the transition region. 
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For the 8^ x 32 zero temperature scale setting calculations with /5=2.0, ci=- 
0.331, mo=1.9, Ls=24 and mj=0.02 58 , the dimensionless lattice meson masses 
are mp=0.95(3) and m^=0. 475(7) and, using nip to set the scale gives a physical 
temperature T=202(6) MeV. Note that mT^/mp=0.5, within errors, for both (3=1.9 
and (3=2.0 simulations. Using (j3.6|) to extrapolate the scale from (3=1.9 gives a 
physical temperature of Text=186 MeV. Violations of asymptotic scaling are around 
eight percent for a five percent change in (3. So, determining the scale from the 
asymptotic scaling relation ()3.6|1 tends to make the transition region appear narrower 
than it would be if a non-perturbative scaling relation were used. 

While there were many interesting indications that observations of true critical 
behavior might be possible on Lt = 4 lattices with domain wall fermions, it is clear 
from the size of the pion mass in physical units that Lg = 24 is too small to suppress 
the large chiral symmetry breaking effects at strong coupling. In section we will 
address the issue of exactly how large an Lg is needed to substantially reduce the mass 
of the dynamical pions in the broken phase and transition region. However, we can 
also tell from our preceding discussion that chiral symmetry breaking effects are much 
smaller in the symmetric phase. Not only does the symmetric phase occur at weaker 
coupling, but reduced average CG iterations per HMC step and the increased rate of 
exponential decay as — » oo for {qq) seen in section 13.21 are strong indications that 
Lg = 24 or even = 16 may be sufficient for simulations in the chirally symmetric 
phase of QCD on = 4 lattices. 



3.4 Effects of finite Ls in the confined phase 

One of the primary advantages of using domain wall fermions to study two flavor QCD 
thermodynamics is the extra parameter Lg which allows us to control chiral symmetry 
breaking. In section 13.11 we concluded that the mass of the pion state indicated that 
the amount of chiral symmetry breaking for domain wall fermions with Lg = 24 was. 
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at best, only moderately better than existing calculations, which are computationally 
less demanding by an order of magnitude. See figure 2 of Karsch's LATTICE 99 
review talk [5^ . 

In defense of the domain wall method, it is important to note that it is currently 
not possible to simulate in the (aTc)~^=4 region of the parameter space with three 
light pions close to their physical mass scale for any of the standard fermion actions. 
So, the relevant questions are why is the pion so heavy for Ls=24 and how large do 
we need to make Ls to reduce all three pions to their physical range. After addressing 
these questions, we will have established a framework to study what improvements are 
possible to reduce chiral symmetry breaking without increasing Lg. The framework 
we propose is the Gell-Mann, Oakes, and Renner relation ()1.85|) discussed in section 
IT31 

In conjunction with the low temperature phase study of the chiral condensate in 
8'^ X 4 volumes with (3=5.2 and mo=1.9, described in section we also measured 
(qq) and x-k at the dynamical mass and various valence masses for several values of 
Ls- Tables fTH and summarize these observables. 

In figure |^ we show the mid-plane contribution to the axial current, AJ5 = 
(qq) — mfX-TT, from ()1.85|) . Since the GMOR relation ()1.85|) holds configuration by 
configuration, we could define a AJ5 for each configuration. However, in our simu- 
lations we only measured the chiral condensate averaged over the spacetime volume 
(qq) and not the contribution from a single point (q^qo) required in ()1.85|) . How- 
ever, the expectation values of two condensates will agree in the ensemble average, 
(qq) ^ {qoqo), our definition of A J5 is valid as the difference of two expectation values. 
We used the jackknife method to estimate the error of AJ5. The long dashed line in 
the figure is the fit AJ5 = 0.0096(2) exp [-0.0190(8)L,] over the range G [16,40] 
with x^/dof=6.5/2. Although the of the fit is a little high, we see that the points 
appear to scatter along the straight line on the log scale. So, our data are consistent 
with chiral symmetry restoration in the Lg 00 limit. 
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We can decompose the effects of finite A J5 by extracting tlie additive mass renor- 
malization m^es and tlie cfiiral condensate subtraction 60 using tlie valence data in 
tables 121] and 121] and the ansatz 

(qq) = {mf + mres) Xn + bo (3.7) 

from p.85j) and p.95|l . We simultaneously fit (qq) and Xn over mj™^"* G [0.02, 0.14] for 
each Lg. To perform this fit, we must further assume Xn can be parameterized as in 
p.88|) . We then fit m^es and Bq to a decreasing exponential in L^. The whole fitting 
procedure is run under a jackknife loop. The fit did not account for correlations 
between data points because there were insufficient data to accurately determine the 
correlation matrix. Our previous experience with this fit ansatz using quenched data 
[22] suggests that including the effects of correlations does not change the answer 
significantly. 

The results are in table |2S1 all errors quoted are jackknife errors. In figure 021 we 
show the L^-dependence of mres and —60 using ten jackknife blocks. The exponential 
fits shown, with jackknife errors, are 

-bo = 0.0106(8) exp [-0.016(3)L,] , xV^of = 0.19(11) (3.8) 
mres = 0.18(1) exp [-0.026(3)L,] , xV^of = 0.08(9) (3.9) 

The fitting range is Lg € [10, 40] and there are four degrees of freedom. Although 
the rate of exponential decay is slow, again we note that our data are consistent 
with chiral symmetry restoration in the Lg —>■ 00 limit. If we use the asymptotic 
scaling relation ()3.6|) to extrapolate the lattice spacing from the 8^ x 32, (3=5.325 
scale setting calculation described in section 13.31 we get ~ 550 MeV with an 
unknown systematic error due to violations of asymptotic scaling. So, in physical 
units, the bare lattice quark mass is m^^^^^'' ^ 10 MeV and the bare residual quark 
mass at L,=24 is m^L*"^'^ ^ 50 MeV. 

In table (211 we show the valence GMOR parameters mres and bo for the scale 
setting calculations |58j which supported the transition region study in section 13.11 
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In the first column, "W" implies the Wilson gauge action ()1.23|) and "R" implies the 
rectangle action (jl.25|) with the choice C\ = —0.331. As before, there were insufficient 
data to accurately estimate the correlation matrix, so it was not included in the fit. 
From the 8^ x 32 scale setting calculations described in section 13. 3t we know that the 
lattice spacing for domain wall fermions with the Wilson gauge action at (3 = 5.325 is 
a~^=652(17) MeV. Using the lattice spacing, we express the dynamical bare lattice 
quark mass as mj''^'^^'*=13.0(3) MeV, in physical units, and the bare residual quark 
mass as mres^'*''^''=37(2) MeV. When compared with the (3=5.2 simulations above, we 
find clear evidence supporting the perturbative expectation that chiral symmetry 
breaking effects are smaller at weaker coupling, for a fixed ^^=24. 

Now, we can estimate whether a fourfold increase in Ls to Ls=9Q would be 
sufficient to reduce the pion masses to a physically interesting scale. First, we will 
imagine that a zero temperature calculation would be done on 8^ x 32 with /5=5.325, 
mo=1.9 and Ls=96. A conservative estimate of the exponential decay of the residual 
mass at /? = 5.325 is mres oc exp(— 0.03Ls) since this agrees with the observed rate 
at (3=5.2 within error. With this decay rate, the bare lattice would be mres ~ 0.007 
at Ls=9Q. Choosing the bare lattice quark mass of m/=0.01 gives m/ > mres and 
will reduce the possibility of large changes in the effective bare mass while scanning 
through the transition region. In [SH], the dynamical bare quark mass dependence 
of the pion mass was determined to be m^ oc 5.38m j when Ls=2A. Assuming the 
same slope will be found when Ls=96 and that the pion mass vanishes (in an infinite 
volume) at mj ^ —0.007, then the dimensionless pion mass will be m^ ^ 0.3. Using 
the lattice spacing found at Ls=24 and m/=0.02, the physical pion mass should be 
j^Cphys) ^ MeV at Ls=96, a physically interesting scale. 

The GMOR relation does provide a useful framework to study the chiral symme- 
try breaking effects of finite domain wall separation. However, one must be careful 
when choosing a parameterization of A J5 that the physical assumptions remain valid 
throughout the region of coupling space of interest. Since our assumption was based 
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on the existence of a pion pole due to spontaneous chiral symmetry breaking, the 
method described here is only applicable in the confined phase. 
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Concluding remarks 



We have established that domain wall fermions preserve the full chiral and flavor 
symmetries of QCD. However, these features come at a substantial cost: one and 
sometimes two orders of magnitude increase in computational effort. So, it is im- 
portant to carefully choose problems where the role of chiral symmetry is crucial to 
determining the correct answer to the problem and similar quality results could not 
be obtained by simply redirecting the increased computational effort into working at 
smaller lattice spacings using less demanding formulations, like staggered or Wilson 
fermions. 

The problem we chose to study is the finite temperature phase transition and 
quark-gluon plasma phase of Nf—2 QCD. This is a particularly challenging problem 
because the finite temperature calculations are usually performed at stronger cou- 
plings than comparable calculations at zero temperature. This is also a particularly 
challenging problem for the domain wall fermion formulation as the suppression of 
chiral symmetry breaking effects at stronger couplings requires larger values of Lg. 

Prom this work it seems clear that domain wall fermions work as advertised. The 
hybrid Monte Carlo (HMC) and conjugate gradient (CG) algorithms worked without 
any difficulty beyond the linear increase in computational cost proportional to Lg. 
The only negative aspect of this work was that the computational resources were not 
four times larger, the amount necessary to simulate on 16^4 lattices with dynamical 
pions close to their physical masses. 
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Despite this limitation, several important physical results were reported here. On 
8^ X 4 volumes, extrapolations to the chiral limit, rrif —>■ and Lg — > oo, on either side 
of the crossover region give a clear evidence of the existence of a finite temperature 
phase transition that restores chiral symmetry in the chiral limit. 

For the first time ever p3j, we performed simulations on 16^ x 4 volumes with 
three dynamical pions whose masses were on the order of 400 MeV. Unfortunately, 
simulations with 400 MeV pions are not sufficiently close to the chiral limit to exhibit 
convincing evidence of the order of the chiral transition. Of course, even reducing the 
mass of the dynamical pions to 200 MeV may not be sufficient to uncover evidence 
of the chiral transition. However, pions are not massless in nature, so the width of 
the crossover region with dynamical pions near their physical value is certainly of 
phenomenological interest. 

Finally, from our study of the GMOR relation, we can determine the residual 
chiral symmetry breaking effects of domain wall fermions at finite Lg. Using the 
GMOR relation, we found these effects were characterized in terms of a residual mass 
rrij-es which decreases exponentially in Lg. In the transition region, for Ls=24 and 
m/=0.02 the effective bare quark mass was dominated by finite Lg effects that were 
too large to allow dynamical pions near the physical pion mass. However, we argued 
quite conservatively that a fourfold increase in Lg with rrif =0.01 should be sufficient 
for the effective bare quark mass to be dominated by in the transition region and 
allow for the masses of the three dynamical pions should be near their physical values. 

We believe the full characterization of the Lt = A transition region will be com- 
pletely accessible to the next generation of supercomputers, even without substantial 
improvement of the domain wall action or algorithms. However, many exciting ques- 
tions about the high temperature phase, like the equation of state can still be 
addressed using the current generation of supercomputers. 
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Table 1: Run parameters for PdrriQ = 1.15) search 
vol: 8^ X 4, mo = 1.15, = 12, m/ = 0.1 

HMC traj. len: ^ x 25, CG stop cond: 10"^ 



/3 start # traj acc. {e~^") plaq. (|W^3|) (qq) 



5.45 





100-800 


0.87 


0.99(1) 


0.470(1) 


0.0168(6) 


0.00276(1) 


5.55 





200-800 


0.87 


0.98(1) 


0.4933(6) 


0.023(1) 


0.002916(6) 


5.65 





300-800 


0.87 


1.00(2) 


0.5218(9) 


0.054(6) 


0.00305(2) 


5.75 


D 


300-800 


0.86 


0.986(9) 


0.5571(7) 


0.196(7) 


0.002875(7) 


5.85 


D 


300-800 


0.85 


1.00(1) 


0.5719(7) 


0.234(3) 


0.002881(3) 


5.95 


D 


200-800 


0.87 


0.99(1) 


0.5857(5) 


0.262(2) 


0.002898(3) 






Table 2: 


Run parameters 


for /3c (mo = 


1.4) search 






vol: 8^ X 4, 


mo = 1.4, 


Ls = 12, 


m/ = 0.1 






HMC traj. len: 


io ^ 25, 


CG stop cond: 10"^ 




P 


start 


# traj 


acc. 




plaq. 


(l^sl) 


{qq) 


5.35 





100-800 


0.87 


1.01(1) 


0.4435(7) 


0.0189(8) 


0.00497(1) 


5.45 





100-800 


0.86 


0.99(1) 


0.4630(7) 


0.0215(8) 


0.00522(1) 


5.55 





300-800 


0.86 


1.00(1) 


0.487(1) 


0.032(3) 


0.00539(3) 


5.65 


D 


400-800 


0.86 


1.00(2) 


0.540(2) 


0.180(6) 


0.00457(4) 


5.75 


D 


300-800 


0.85 


0.99(1) 


0.5598(8) 


0.224(3) 


0.00445(1) 


5.85 


D 


200-800 


0.89 


1.011(8) 


0.5744(3) 


0.254(3) 


0.004409(5) 
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Table 3: Run parameters for PdrriQ = 1.65) search 
vol: 8^ X 4, mo = 1.65, Ls = 12, m/ = 0.1 
HMC traj. len: ^ x 25, CG stop cond: 10"^ 



/3 start # traj acc. {e~^") plaq. (IW3I) (qq) 



5.25 





200-800 


0.82 


0.97(2) 


0.4289(5) 


0.027(2) 


0.01000(2) 


5.35 





400-800 


0.68 


0.98(4) 


0.451(3) 


0.035(4) 


0.01000(9) 


5.45 


D 


400-800 


0.74 


1.09(5) 


0.4769(8) 


0.049(7) 


0.00985(7) 


5.55 


D 


600-1200 


0.80 


1.00(3) 


0.531(1) 


0.175(7) 


0.00718(7) 


5.65 


D 


400-800 


0.79 


0.96(2) 


0.5507(9) 


0.216(4) 


0.00677(2) 


5.75 


D 


200-800 


0.88 


0.989(7) 


0.5663(4) 


0.243(3) 


0.00658(1) 



Table 4: Run parameters for PdruQ = 1.8) search 





vol: 8^ X 4, 


mo = 


= 1.8, 


Ls = 12, 


m/ = 


0.1 






HMC traj. 


len: ^ 


X 25, 


CG stop cond: 10 


-6 




start 


# traj 


acc. 




plaq. 


{m) 




5.15 





200-800 


0.83 


1.02(2) 


0.4191(8) 


0.029(1) 


0.01485(5) 


5.25 





400-800 


0.66 


0.97(5) 


0.4381(6) 


0.038(2) 


0.01458(5) 


5.35 





400-800 


0.63 


0.97(5) 


0.471(2) 


0.052(3) 


0.0134(2) 


5.45 





400-800 


0.76 


1.01(3) 


0.515(2) 


0.161(4) 


0.0097(1) 


5.55 


D 


400-800 


0.79 


1.05(5) 


0.540(1) 


0.200(9) 


0.0088(1) 


5.65 


D 


200-800 


0.89 


1.01(2) 


0.5570(5) 


0.242(4) 


0.00828(2) 
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Table 5: Run parameters for (5c{mQ — 1.9) search 







vol: 8^ X 4, mo 


= 1.9, 




: 12, ruf 


= 0.1 










CG stop cond: 10 


6 








start 


traj. len. 


# traj 


acc. 




plaq. 


(^31) 


{qq) 


5.0 





i)X20 


200-800 


0.37 


0.8(1) 


0.4002(8) 


0.032(2) 


0.01919(5) 


5.2 





ix20 


400-800 


0.43 


1.2(3) 


0.437(1) 


0.032(2) 


0.01838(4) 


5.25 





^x25 


400-800 


0.65 


1.10(9) 


0.452(1) 


0.049(6) 


0.0174(2) 


5.35 


D 


^x25 


600-1200 


0.69 


0.95(5) 


0.493(2) 


0.107(9) 


0.0135(4) 


5.45 


D 


^x25 


600-1200 


0.74 


1.01(4) 


0.528(1) 


0.197(4) 


0.01039(7) 


5.55 


D 


MX 25 


400-830 


0.82 


1.00(1) 


0.5463(5) 


0.227(6) 


0.00974(4) 


5.65 


D 


^x25 


400-800 


0.88 


1.03(1) 


0.5613(8) 


0.248(5) 


0.00943(4) 



Table 6: Run parameters for (5c{mQ — 2.0) search 





vol: 8^ X 4, 


mo = 


= 2.0, 


Ls = 12, 


m/ = 


0.1 






HMC traj. 


len: ^ 


X 25, 


CG stop cond: 10"^ 






start 


# traj 


acc. 




plaq. 


m^\) 


{qq) 


5.05 





200-800 


0.77 


0.99(3) 


0.4192(8) 


0.035(1) 


0.02324(8) 


5.15 





200-800 


0.75 


0.98(3) 


0.442(1) 


0.042(3) 


0.0215(2) 


5.25 





200-1200 


0.79 


1.03(1) 


0.474(1) 


0.080(7) 


0.0181(3) 


5.35 


D 


200-800 


0.83 


1.00(2) 


0.5130(7) 


0.173(6) 


0.0130(2) 


5.45 


D 


200-800 


0.87 


1.02(2) 


0.5349(5) 


0.203(3) 


0.01157(4) 


5.55 


D 


200-800 


0.85 


1.01(1) 


0.5503(4) 


0.235(3) 


0.01099(2) 
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Table 7: Run parameters for (5c{mQ = 2.15) search 





vol: 8^ X 4, 


mo = 


-- 2.15, 


Ls = 12, 




0.1 






HMC traj. 


len: ^ 


X 25, 


CG stop cond: lO'^ 




3 


sta.rt 


# tra,i 


ace. 




plaq. 


(m) 


(w) 


4.85 





200-800 


0.69 


0.95(3) 


0.4004(8) 


0.034(2) 


0.0323(2) 


4.95 





200-800 


0.72 


0.97(3) 


0.419(2) 


0.040(2) 


0.0302(3) 


5.05 





200-800 


0.48 


0.92(4) 


0.443(2) 


0.052(3) 


0.0272(5) 


5.15 





200-1200 


0.62 


1.01(3) 


0.480(3) 


0.12(1) 


0.0203(7) 


5.25 





400-800 


0.70 


0.97(5) 


0.5105(4) 


0.185(2) 


0.01559(8) 


5.35 


D 


400-800 


0.69 


1.01(4) 


0.529(1) 


0.216(6) 


0.0141(1) 


5.45 


D 


400-800 


0.71 


1.00(3) 


0.5453(7) 


0.230(4) 


0.01330(5) 



Table 8: Run parameters for PdrriQ — 2.4) search 
vol: 8^ X 4, mo = 2.4, = 12, m/ = 0.1 

HMC traj. len: ^ x 25, CG stop cond: lO'^ 



13 start # traj acc. {e~^^) plaq. (|H^3|) {qq) 



4.65 





100-800 


0.63 


1.02(5) 


0.3953(6) 


0.046(3) 


0.0484(3) 


4.75 





200-800 


0.68 


0.99(5) 


0.4156(9) 


0.054(3) 


0.0442(3) 


4.85 





300-800 


0.70 


0.94(4) 


0.439(2) 


0.069(5) 


0.0380(6) 


4.95 





200-800 


0.77 


1.02(4) 


0.4779(6) 


0.155(4) 


0.0257(2) 


5.05 


D 


200-800 


0.80 


1.01(2) 


0.4987(9) 


0.190(2) 


0.0220(2) 


5.15 


D 


200-800 


0.84 


1.01(3) 


0.5170(5) 


0.221(3) 


0.01962(7) 



Table 9: Fit of (|M^3|) = Cp {ci + tanh [c2((3 - (3c)]}- 8^ x 4, = 12, ruf = 



mo 


/3c 


Co 


Cl 


C2 






1.15 


5.722(5) 


0.122(1) 


1.126(7) 


9.8(5) 


2 


14.50 


1.4 


5.627(4) 


0.110(1) 


1.180(7) 


17(1) 


2 


21.39 


1.65 


5.524(6) 


0.107(2) 


1.24(2) 


10(1) 


2 


9.60 


1.8 


5.431(4) 


0.105(2) 


1.29(1) 


11.7(8) 


2 


11.07 


1.9 


5.386(6) 


0.107(3) 


1.27(2) 


9.9(8) 


3 


2.61 


2.0 


5.327(7) 


0.107(3) 


1.27(3) 


6.7(6) 


2 


10.33 


2.15 


5.179(6) 


0.098(2) 


1.35(3) 


8.3(6) 


3 


1.13 


2.4 


4.934(7) 


0.091(4) 


1.46(5) 


7.2(6) 


2 


14.02 



10: 


Fit of (qq) 


= Co {ci + tanh [c2(/3 - 


Pc)]}: 8' 


X 4, Ls 


= 12, m/ 


mo 




Co 


Cl 


C2 


Ndoi 




1.65 


5.524(6) 


-0.00170(1) 


-4.89(3) 


18(2) 


2 


38.28 


1.8 


5.398(4) 


-0.00331(5) 


-3.50(3) 


10.6(7) 


2 


4.90 


1.9 


5.330(3) 


-0.00488(3) 


-2.94(2) 


9.1(2) 


3 


2.85 


2.0 


5.256(5) 


-0.00648(9) 


-2.67(3) 


7.5(3) 


2 


7.12 


2.15 


5.103(6) 


-0.0099(1) 


-2.32(2) 


6.5(3) 


3 


3.36 


2.4 


4.864(4) 


-0.0154(3) 


-2.24(3) 


6.9(3) 


2 


14.96 
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Table 11: Run parameters for 8^ x 4,/? = 5.2, mo = 1.9 confined phase study 



Lg traj. len. # traj. 



acc. (e ^■'^) plaq. 



n no A 
yj.yjz 4 


1 

64 


X OZ 


ZUU-oUU 


n Q/i 
u.y4 


1 nno/'a^ 

i.UUZ(^o ) 






yj.yjzoyjy i ) 




1 

64 


X 32 


200-800 

^yjyj Kjyjyj 


Q2 

yj. u ^ 


1 OOfl 

-L . \J\J \ ^ ) 


468('2'l 


QRQ(A) 

W . WOW I t: I 


0167('3'l 

W . W -L W 1 1 


8 

o 


1 

64 


X 32 


400-800 


8Q 


QSfl 


456l'2'l 


061 (2) 

W . W W X 1 ^ ( 


01 33('2'1 


1 n 


1 

64 


X 32 


200-2000 


8fi 


QQ5('8'l 




048l'2'l 

W. WtiO I Z( J 


01 24l'l 


12 


1 

64 


X 32 


200-2000 

^yjyj ^yjyjyj 


84 


1 01 fi 1 


n 4428f6l 


048f2l 


01 123f7l 


Ifi 


1 

64 


X 32 


550-2000 


75 


Q8/'2'l 


4388fQ'l 




00Q87fQ'l 


94 


1 

100 


X "lO 


350-2000 


73 


Q5('2'l 


435Q('7'l 


047('3'l 


0088^1 


32 


1 

100 


X 50 


300-2000 


72 


1 03('2') 


\'W1(1\ 

\) . J. 1 1 i 1 


045('2') 


00835 ('7'! 

W.WWOtJtJl 1 I 


4(1 


1 

128 


X 64 


300-1 350 


73 


1 0?('3'l 




n 044('2'l 


Ci077?(R) 


u.uu o 


1 

50 


X ZD 


zuu-you 


U.oO 


u.yo^^i J 






U.Ui ( uyz J 


16 


1 

64 


X 32 


200-820 


0.84 


0.99(2) 


0.4361(8) 


0.045(3) 


0.0135(1) 


0.1 4 


1 

40 


X 20 


300-800 


0.82 


1.00(2) 


0.463(2) 


0.058(3) 


0.0375(4) 


6 


1 

50 


X 25 


300-800 


0.89 


1.00(2) 


0.454(3) 


0.048(5) 


0.0246(4) 


8 


1 

40 


X 20 


300-800 


0.57 


0.89(5) 


0.4437(8) 


0.040(3) 


0.02109(7) 


10 


1 

50 


X 25 


200-800 


0.83 


1.00(3) 


0.4405(6) 


0.036(2) 


0.01927(9) 


12 


1 

40 


X 20 


400-800 


0.43 


1.2(3) 


0.437(1) 


0.032(2) 


0.01838(4) 


16 


1 

64 


X 32 


200-800 


0.80 


0.98(2) 


0.435(2) 


0.035(2) 


0.01709(9) 


24 


1 

64 


X 32 


200-800 


0.72 


0.94(2) 


0.433(1) 


0.033(2) 


0.01596(7) 


32 


1 

100 


X 50 


200-800 


0.82 


0.99(2) 


0.4305(5) 


0.037(2) 


0.01547(7) 


40 


1 

100 


X 50 


200-800 


0.78 


1.00(5) 


0.432(1) 


0.035(2) 


0.01524(5) 



0.14 8 

16 



i)X20 

^x32 



200-860 

200-800 



0.63 1.03(6 
0.85 1.03(2 



0.4433(7) 

0.433(1) 



0.033(5 
0.030(1 



0.0241(1) 

0.02017(6) 



0.18 8 
16 



^x20 
^x32 



200-1200 
200-800 



0.70 1.02(3 
0.84 0.98(2 



0.4410(7) 
0.432(1) 



0.030(1 
0.033(1 



0.02686(7) 
0.02309(5) 



Table 12: 


Fit of (qq) = Cq + Ciuif + C2m^, 


8^ X 4,/5 


= 5.2, mo = 1.9 




fit range 


Co 


Cl 


C2 


A^dof 




8 


0.02- 


-0.1 


0.0117(2) 


0.095(2) 




1 


3.67 




0.02- 


-0.14 


0.0122(2) 


0.088(2) 




2 


9.89 




0.02- 


-0.18 


0.0130(1) 


0.0781(9) 




3 


20.48 




0.02- 


-0.14 


0.0111(3) 


0.118(7) 


-0.18(4) 


1 


0.23 






-0 1 R 

U. lO 


0.0112(3)'^ 


0.114(5) 


-0.15(2) 


2 


0.37 


10 


0.02- 


-0.1 


0.0107(2)" 


0.086(2) 









12 


0.02- 


-0.1 


0.00944(9)^^ 


0.089(1) 









16 


0.02- 


-0.1 


0.0081(1) 


0.090(2) 




1 


0.24 




0.02- 


-0.14 


0.00830(9) 


0.0855(9) 




2 


6.68 




0.02- 


-0.18 


0.00857(8) 


0.0815(6) 




3 


16.38 




0.02- 


-0.14 


0.0079(2) 


0.101(5) 


-0.10(3) 


1 


0.73 




0.02- 


-0.18 


0.0079(1)" 


0.100(3) 


-0.09(1) 


2 


0.43 


24 


0.02- 


-0.1 


0.0069(1)" 


0.090(2) 









32 


0.02- 


-0.1 


0.00656(9)" 


0.089(1) 









40 


0.02- 


-0.1 


0.00584(9)" 


0.094(1) 











0.02- 


-0.1 


0.0060(1) 


0.092(1) 










"Shown in figure I?!! and used as input to fits in table 1131 
''Used fitted values from table 1131 as input. 
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Table 13: 


Fit of {qq) = Cq + Ci exp 


) 8^ X 4 8 


= 5.2, 


777,Q =1.9 


rrif 


fit 


range 


Co 


Cl 




' Qui 


A / ClOl 




8- 


-40 


0.0059(1) 


0.014(1) 


0.11(1) 


4 


6.62 




10- 


-40 


0.0060(1) 


0.016(2) 


0.13(1) 


3 


7.90 




12- 


-40 


0.0058(2) 


0.012(2) 


0.10(1) 


2 


9.41 




16- 


-40 


0.003(4) 


0.007(3) 


0.02(2) 


1 


3.16 


0.02 


8- 


-40 


0.00779(8)^ 


0.014(1) 


0.116(8) 


4 


5.55 




10- 


-40 


0.00780(9) 


0.014(1) 


0.118(9) 


3 


7.37 




12- 


-40 


0.0076(1) 


0.011(1) 


0.10(1) 


2 


7.74 




16- 


-40 


0.0065(9) 


0.0064(4) 


0.04(2) 


1 


4.67 


0.1 


8- 


-40 


0.01527(4) 


0.0188(8) 


0.149(5) 


4 


5.14 




10- 


-40 


0.01514(6)^ 


0.013(1) 


0.119(8) 


3 


0.41 




12- 


-40 


0.01513(7) 


0.013(1) 


0.117(9) 


2 


0.56 




16- 


-40 


0.0151(1) 


0.010(3) 


0.101(2) 


1 


0.03 



"Used fitted values from table 1121 as input. 
''Used as input to fits in tablell2l 
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Table 14: Run parameters for 8^ x 4, f3 = 5.45, mo = 1.9 deconfined phase study 





Ls 


traj. len. 


# traj. 


acc. 




plaq. 




im) 


0.02 


8 


1 

64 


X 32 


200-800 


0.91 


1.005(7) 


0.5376(7) 


0.226(4) 


0.00415(6) 




10 


1 

64 


X 32 


200-1000 


0.91 


0.992(8) 


0.5328(6) 


0.207(4) 


0.00319(5) 




12 


1 

64 


X 32 


200-800 


0.95 


1.009(9) 


0.5300(4) 


0.202(5) 


0.00270(3) 




16 


1 

64 


X 32 


200-800 


0.90 


1.02(1) 


0.5266(8) 


0.199(4) 


0.00237(6) 




24 


1 

64 


X 32 


400-1200 


0.86 


0.98(2) 


0.5257(7) 


0.187(3) 


0.00216(6) 




32 


1 
100 


X 50 


400-800 


0.94 


1.00(2) 


0.524(2) 


0.180(5) 


0.00209(5) 


0.06 


8 


1 

50 


X 25 


200-1000 


0.86 


0.99(3) 


0.536(1) 


0.217(3) 


0.0080(1) 




10 


1 

64 


X 32 


200-1000 


0.92 


0.994(7) 


0.5313(6) 


0.203(4) 


0.00704(5) 




12 


1 

64 


X 32 


200-1000 


0.89 


1.013(8) 


0.5286(8) 


0.195(4) 


0.00666(5) 




16 


1 

64 


X 32 


400-800 


0.76 


1.02(4) 


0.525(2) 


0.192(4) 


0.00637(7) 




24 


1 

64 


X 32 


300-1000 


0.84 


1.00(1) 


0.521(2) 


0.174(6) 


0.00617(9) 




32 


1 

64 


X 32 


500-1000 


0.80 


1.00(2) 


0.525(2) 


0.189(3) 


0.00592(4) 


0.1 


8 


1 

50 


X 25 


300-800 


0.83 


0.98(2) 


0.5336(6) 


0.211(4) 


0.01174(4) 




10 


1 

50 


X 25 


300-990 


0.88 


0.99(1) 


0.5310(9) 


0.200(2) 


0.01075(5) 




12 


1 

50 


X 25 


600-1200 


0.74 


1.01(4) 


0.528(1) 


0.197(4) 


0.01838(4) 




16 


1 

64 


X 32 


400-800 


0.79 


1.01(3) 


0.523(1) 


0.170(5) 


0.0103(1) 




24 


1 

64 


X 32 


400-2000 


0.86 


0.991(8) 


0.512(1) 


0.170(8) 


0.0102(1) 




32 


1 

()i 


X 32 


300-1000 


0.81 


0.98(2) 


0.519(1) 


0.159(5) 


0.01011(9) 


0.14 


8 


1 

50 


X 25 


200-800 


0.83 


1.01(1) 


0.533(1) 


0.210(3) 


0.01531(9) 




16 


1 

64 


X 32 


600-1200 


0.76 


0.98(2) 


0.520(1) 


0.159(9) 


0.0143(1) 


0.18 


8 


1 

50 


X 25 


400-800 


0.81 


1.03(2) 


0.5314(6) 


0.202(4) 


0.01884(5) 




16 


1 

64 


X 32 


600-1200 


0.78 


0.94(2) 


0.515(1) 


0.141(8) 


0.0182(2) 
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Table 15: Fit of {qq) = Cq + Ciuif + C2m^, 8^ x 4, /3 = 5.45, mo = 1.9 





fit range 


Co 


Cl 


C2 


A^dof 


xViVdof 


8 


0.02-0.1 


0.00227(7) 


0.0948(9) 


— 


1 


0.62 


0.02-0.18 


0.00219(9)" 


0.099(2) 


-0.037(9) 


2 


0.10 


0.00246(6) 


0.0915(5) 


— 


3 


5.80 


10 


0.02-0.1 


0.00131(6)" 


0.0947(9) 


— 


1 


1.23 


12 


0.02-0.1 


0.00077(4)" 


0.0968(8) 


— 


1 


3.85 


16 


0.02-0.1 


0.00039(8) 


0.100(2) 





1 


0.09 





0.1060(8) 





2 


12.86 


0.02-0.18 


0.0004(1) 


0.100(3) 


-0.01(2) 


2 


0.02 


0.00040(6)" 


0.0992(8) 


— 


3 


0.07 




0.111(1) 


-0.06(1) 


3 


4.81 


24 


0.02-0.1 


0.00016(8)" 


0.100(2) 




1 


0.01 




0.103(1) 


-0.07(3) 


1 


0.83 




0.103(1) 




2 


2.20 


32 


0.02-0.1 




0.100(5) 




2 


3.99 


0.00006(6)" 


0.099(1) 




1 


7.11 




0.099(2) 


0.02(2) 


1 


7.40 




0.02-0.1 




0.1012(4) 




2 


6.90 


0.00011(5) 


0.0995(8) 




1 


8.45 



"Shown in figure 1^ and used as input to fits in table 1161 
''Used fitted values from table 1161 as input. 
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Table 16: Fit of {gg) = cp + Ci exp (-caL^), 8^ x 4, /? = 5.45, mp = 1.9 





fit range 


Co 


Cl 


C2 


A^dof 




0.0'^ 


8-32 


0.00010(5) 


0.019(3) 


0.27(2) 


3 


0.76 








0.015(2) 


0.24(1) 


4 


1.64 




10-32 


0.00013(5) 


0.02(1) 


0.30(4) 


2 


1.45 








0.011(3) 


0.22(2) 


3 


1.74 




12-32 




0.005(2) 


0.15(3) 


2 


0.22 






0.00014(5) 


0.03(5) 


0.3(1) 


1 


4.32 


0.02 


8-32 


0.00213(4)^ 


0.025(4) 


0.31(2) 


3 


1.07 




10-32 


0.00212(4) 


0.018(6) 


0.28(4) 


2 


1.24 




12-32 


0.00208(6) 


0.005(3) 


0.18(6) 


1 


0.10 


0.06 


8-32 


0.00599(4)^ 


0.015(3) 


0.26(2) 


3 


4.84 




10-32 


0.00591(6) 


0.005(1) 


0.15(3) 


2 


2.05 




12-32 


0.00598(5) 


0.01(1) 


0.24(7) 


1 


8.04 


0.1 


8-32 


0.01016(6)^ 


0.08(3) 


0.48(6) 


3 


0.42 




10-32 


0.01015(7) 


0.04(7) 


0.4(2) 


2 


0.58 




12-32 


0.0101(1) 


0.002(6) 


0.2(3) 


1 


0.30 



"Used fitted values from table 1151 as input. 
''Used as input to fits in tablellSI 
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Table 17: Run parameters for transition study: 8^ x 4, mo = 1.9, = 24, m/ = 0.02 

13 start traj. len. # traj. acc. (e.~^^) plaq. (iW^sl) 

5.2 O 1^x50 300-2000 0.72 0.95(2) 0.4359(6) 0.046(2) 0.0087(1) 

5.275 O T5o^50 200-1000 0.83 1.01(2) 0.459(1) 0.068(4) 0.0069(2) 

5.325 D 1^x50 200-1150 0.89 1.01(1) 0.480(3) 0.10(1) 0.0050(4) 

5.45 D ^x32 400-1200 0.86 0.98(2) 0.5257(7) 0.187(3) 0.00216(6) 



Table 18: Run parameters for transition study: 16'^ x 4, mg = 1.9, Lg = 24, mf = 0.02 







HMC traj. len: ^ 


X 64, 


CG stop cond: 10"^ 




/? 


sta.rt 


#t 


raj. 


acc. 


(c---> 


plaq. 




(w) 


5.25 


D 


140- 


-440 


0.69 


1.01(6) 


0.4465(2) 


0.046(2) 


0.0087(1) 







350- 


-650 


0.74 


0.99(5) 


0.4475(3) 


0.051(2) 


0.00814(4) 


5.275 


D 


360- 


-660 


0.70 


0.92(5) 


0.4565(5) 


0.062(2) 


0.00737(4) 







470- 


-770 


0.71 


0.93(5) 


0.4561(3) 


0.061(2) 


0.00747(4) 


5.3 


D 


1100- 


-1400 


0.72 


1.05(4) 


0.4656(8) 


0.074(2) 


0.0066(1) 







720- 


-1020 


0.75 


1.06(4) 


0.4683(6) 


0.077(2) 


0.00629(7) 


5.325 


D 


490- 


-790 


0.82 


1.00(3) 


0.4800(6) 


0.096(2) 


0.00506(6) 







670- 


-970 


0.77 


0.96(3) 


0.480(1) 


0.102(4) 


0.0049(3) 


5.35 


D 


650- 


-950 


0.86 


1.04(3) 


0.4924(9) 


0.124(3) 


0.0037(1) 







200- 


-500 


0.84 


1.03(3) 


0.4928(4) 


0.127(2) 


0.00366(5) 
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Table 19: Run parameters for 8^ x 4 improved gauge action transition study 



vol: 


8=^ X 4, 


Cl = 


-0.331 


mo 


= 1.9, 


Ls = 24, 


ruf = 0.02 






HMC traj. 


len: ^ x 50, 


CG stop cond: 10 


-6 




start 


^ traj. 


acc. 




plaq. 


(m) 


\ 1 1 / 


1.7 


D 


70-170 


0.77 


1.00(5) 


0.424(1) 


0.032(2) 


0.01067(9) 







70-170 


0.76 


0.96(4) 


0.423(1) 


0.034(2) 


0.0104(2) 


1.8 


D 


100-200 


0.86 


1.07(5) 


0.4620(9) 


0.037(2) 


0.0097(2) 







100-200 


0.78 


0.93(4) 


0.462(1) 


0.033(4) 


0.0096(2) 


1.9 


D 


130-230 


0.88 


0.98(5) 


0.507(2) 


0.057(3) 


0.0066(2) 







130-230 


0.89 


1.11(6) 


0.5075(9) 


0.064(5) 


0.0064(2) 


2.0 


D 


160-260 


0.93 


0.98(3) 


0.5500(8) 


0.123(6) 


0.0028(1) 







160-260 


0.96 


0.99(3) 


0.5502(9) 


0.128(3) 


0.00245(5) 


2.1 


D 


200-400 


0.94 


0.98(1) 


0.5824(6) 


0.17(1) 


0.0021(2) 







200-400 


0.94 


0.99(3) 


0.5806(7) 


0.17(1) 


0.0023(2) 


2.2 


D 


140-240 


0.95 


1.01(1) 


0.6054(4) 


0.20(2) 


0.00191(3) 







140-240 


0.93 


0.99(1) 


0.6064(5) 


0.212(9) 


0.00187(1) 
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Table 20: Run parameters for 16^ x 4 improved gauge action transition study 
vol: 163 X 4, ci = -0.331, mo = 1.9, = 24, m/ = 0.02 

HMC traj. len: ^ x 64, CG stop cond: 10"^ 
P start # traj. acc. (e"^^) plaq. (|W^3|> (qq) 



1.85 


D 


150- 


-270 


0.75 


1.0(1) 


0.4848(4) 


0.044(1) 


0.00843(6) 







160- 


-300 


0.73 


1.09(9) 


0.4861(4) 


0.044(2) 


0.00823(7) 


1.9 


D 


200- 


-300 


0.71 


0.96(6) 


0.5079(7) 


0.063(2) 


0.0064(1) 







200- 


-300 


0.78 


0.99(3) 


0.5095(3) 


0.0658(8) 


0.00598(7) 


1.95 


D 


220- 


-520 


0.76 


1.00(2) 


0.5313(3) 


0.098(2) 


0.00378(5) 







150- 


-450 


0.83 


1.04(4) 


0.5321(3) 


0.107(2) 


0.00362(4) 


2.0 


D 


200- 


-300 


0.93 


1.08(3) 


0.5506(3) 


0.1242(6) 


0.00260(7) 







200- 


-300 


0.90 


1.05(3) 


0.5505(3) 


0.127(2) 


0.00260(5) 
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Table 21: Valence {qq) and x^- 8^ x 4, /? = 5.2, mp ^ 1.9, = 10-16, m/ = 0.02 







{m) 




r (dyn) 


(val) 


# traj 


# meas 


avg 


# traj 


# meas 


avg 


10 


0.02 


200-2000 


1800 


0.0124(1) 


800-2000 


240 


0.130(1) 




0.06 




360 


0.0156(1) 






0.1202(7) 




0.1 






0.0186(1) 






0.1130(6) 




0.14 






0.0215(1) 






0.1071(5) 




0.18 

0.22 






0.02432(9) 

0.02695(8) 






0.1022(4) 

0.0979(4) 


12 


0.02 


200-2000 


1800 


0.01123(7) 


800-2000 


240 


0.129(1) 




0.06 




360 


0.01448(8) 






0.1193(7) 




0.1 






0.01759(7) 






0.1120(6) 




0.14 






0.02056(6) 






0.1061(4) 




0.18 






0.02338(5) 






0.1012(4) 




0.22 






0.02608(5) 






0.0970(3) 


16 


0.02 


550-2000 


1450 


0.00987(9) 


800-2000 


240 


0.134(1) 




0.06 




290 


0.01318(8) 






0.123(1) 




0.1 






0.01639(7) 






0.1146(8) 




0.14 






0.01945(6) 






0.1082(7) 




0.18 






0.02236(6) 






0.1029(6) 




0.22 






0.02513(5) 






0.0985(5) 
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Table 22: Valence {qq) and Xn- 8^ x 4, /3 = 5.2, mp = 1.9, L, = 24-40, m/ = 0.02 







{qq) 




^((lyn) 


(val) 

III 

J 


4i trni 
it ' -"-f'.l 


-jj- on. s 




^ trni 
it ^^".\ 






0/1 


n no 


obU— zUUU 


iboU 




onn onnn 
oUU— zUUU 


o/i n 
z4U 


n 1 /I T /^o^ 
U.141(^Z j 




u.uo 




330 


0.0122(1) 






0.127(1) 




0.1 






0.01550(9) 






0.1177(9) 




0.14 






n ni fiRQ/'s^ 

U.UioDol^o j 






nil ns/'7^ 




U. iO 






n noi ^1 
U.Uzioil^ / ) 






n int^o/'t^A 




n 00 






n no/i /I 
U.UZ44ol^0 j 






n 1 nnc;^/i'> 
U.lUU0t^4 ) 




n no 


oUU— zUUU 




n nnQQ /" 
U.UUooo( / ) 


Qnn onnn 
oUU— zUUU 


o/i n 
z4U 


U.ioo(o ) 




n nfi 

U.Uu 




340 


0.01169(8) 






0.135(2) 




0.1 






0.01504(8) 






0.123(2) 




0.14 






0.01822(7) 






0.114(1) 




0.18 






0.02124(6) 






0.108(1) 




0.22 






0.02410(6) 






0.103(1) 


40 


0.02 


300-1350 


1050 


0.00772(8) 


300-1350 


210 


0.157(3) 




0.06 




210 


0.0113(1) 






0.135(2) 




0.1 






0.0147(1) 






0.123(1) 




0.14 






0.01799(7) 






0.114(1) 




0.18 






0.02106(6) 






0.1075(9) 




0.22 






0.02396(6) 






0.1022(7) 
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Table 23: Valence m^es and 60: 8^ x 4,/3 = 5.2, mp = 1.9, m/ = 0.02 













r 


blks 


' ' 'res 


I'D 


/dof 


1 n 


1 n 
iU 


u.i4y(^o j 


n nno/1 (k\ 
U.uuy4(^0 ) 


U. ( 0(^4Z J 




io 


u. i4y (^0 J 


n nnQ/i 
u.uuy4(^o ) 


U. ( ol^4i J 




on 


n 1 A(\( A \ 


U.uuy4\^4 ) 


n 7Q/'/i /I ^ 
U. / o\^44 J 


1 


1 n 


n 1 oo/'o^ 


n nnsn/'o^ 

U.UUoUl^Z ) 


1 ^^Q/'/l/l^ 
i.Do(^44 J 




io 


u. izy (^z J 


n nnfin('/i ^ 
U.UUoU(^4 ) 


i.0o(^0y j 




on 


n 1 oo^o^ 
u.izy(^z j 


n nnsn^Q^ 


1 7n^7n^ 
1. lUy / U J 


ID 


1 n 
iU 




U.UUoU\^4 ) 


1 1 /I /'/I 
1.14(^40 J 




io 


nil Q^'Q'i 
u. i io(^o j 


n nnsn^'/i 

U.UUoU(^4 j 


1 on^''?7'\ 




on 


nil "^l A \ 

U.iiOl^4 ) 


n nnsn^'c^'i 

U.UUoUl^O J 


i.Zi (^4D J 


24 


10 


0.095(2) 


0.0075(3) 


1.52(74) 




15 


0.095(3) 


0.0075(3) 


1.57(68) 




20 


0.095(3) 


0.0075(3) 


1.61(73) 


32 


10 


0.078(2) 


0.0068(5) 


0.71(38) 




15 


0.078(2) 


0.0069(4) 


0.72(39) 




20 


0.078(2) 


0.0069(4) 


0.72(42) 


40 


10 


0.059(3) 


0.0048(3) 


1.70(93) 




15 


0.059(2) 


0.0048(3) 


1.74(1.47) 




20 


0.060(2) 


0.0048(3) 


1.74(1.56) 



100 



Table 24: Valence m^cs and bo for scale setting simulations, tuq = 1.9 



gauge 
action 


Vol 




Ls 


ruf 


jknf 
blks 




-bo 


xVdof 


W 


8x32 


5.325 


24 


0.02 


10 


0.057(2) 


0.0045(2) 


8.4(3.9) 


20 


0.057(2) 


0.0046(2) 


8.8(3.9) 


0.06 


10 


0.059(3) 


0.0045(3) 


2.1(6) 


20 


0.059(2) 


0.0045(2) 


2.2(1.0) 


R 


8x32 


1.9 


24 


0.02 


10 


0.038(2) 


0.0027(1) 


8.6(2.2) 


20 


0.038(2) 


0.0027(2) 


9.1(2.8) 


2.0 


24 


0.02 


10 


0.016(2) 


0.0014(2) 


6.8(2.6) 


20 


0.016(3) 


0.0015(4) 


7.0(5.1) 


0.06 


10 


0.017(2) 


0.0015(3) 


1.6(8) 


20 


0.017(2) 


0.0015(3) 


1.6(8) 
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Figures 
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Figure 2: (11^31) evol: 8^ x 4,mo = 1.15, = 12, mj = 0.1 
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Figure 4: (11^31) evol: 8^ x 4, mo = 1.4, = 12, m/ = 0.1 
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Figure 5: (gg) evol: 8^ x 4, mo = 1.65, = 12, m/ = 0.1 
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Figure 6: (11^31) evol: 8^ x 4,mo = 1.65, = 12, m/ = 0.1 
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Figure 8: {{WsD evol: 8^ x 4, mo = 1.8, = 12, m/ = 0.1 
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Figure 9: (qq) evol: 8^ x 4, mo = 1.9, = 12, m/ = 0.1 
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Figure 11: {qq) evol: 8^ x 4, mo = 2.0, = 12, m/ = 0.1 
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Figure 13: {qq) evol: 8^ x 4, mo = 2.15, = 12, m/ = 0.1 
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Figure 14: {{WsD evol: 8^ x 4, mo = 2.15, = 12, mj = 0.1 
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Figure 15: {qq) evol: 8^ x 4, mo = 2.4, = 12, m/ = 0.1 
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Figure 16: {{WsD evol: 8^ x 4, mo = 2.4, = 12, m/ = 0.1 
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Figure 17: tanh fits of I3c. 8^ x 4; mo = 1.15, ■ ■ ■ , 1.8; = 12; m/ = 0.1 
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Figure 19: f3c versus mo: 8^ x 4, = 12, mf = 0.1 
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Figure 20: Fit of {qq) = Cq + CiTUf + CsmJ: 8^ x 4; (3 = 5.2; mo = 1.9; L,, = 8, 16 
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Figure 21: Fit of (qq) = Cq + Ci exp (-Ca-L,): 8^ x 4, /3 = 5.2, mo = 1.9 
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Figure 22: Fit of {qq) = Cq + Cirrij + CsmJ: 8^ x 4; /5 = 5.45; ttiq = 1.9; = 8, 16 



124 



♦ fn,=0.1 
A m,=0.06 
i m,=0.02 

0.02 - 1 V m,^0 



0.005 



1 

1 

\ 

la 

V 



\ 

0.01 - \ 



\ 



Q — I — I — I — \ — I — I — I — 'v " 1 — I — I — — 1 — I — L \ \ \ 

8 16 24 32 40 

L. 



Figure 23: Fit of (qq) = co + ci exp (-C2-Ls): 8^ x 4, /3 = 5.45, mo = 1.9 
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Figure 26: (gg), {{Wsl) vs. (3: 8^ x 4, mo = 1.9, = 24,m/ = 0.02 
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Figure 27: (qq) evol: 16^ x 4, mo = 1.9, = 24, m/ = 0.02 



0.1 - 



V ' I' I 

H I 



I I I 



I I I 



I I I I I I I 



II ill 11 




J I b I I I I I I I I I I I I I I I I I I I I I I I L 



0.1 - 







-r — I — I — I — I — I — q — I — n— 1 — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — r 

. .. I I 



p=5.275 . 



11 
I I I I I I I I I I I I I I I I I L 



0.1 - 



TTT 1 1 1 1 1 1 1 1 1 1 1 TT— I 1 1 1 1 1 1 1 1 1 1 1 r 




0.1 - 



I'T^I 1 1 1 1 1 1 1 


ll 1 1 ll 1 


1 1 1 1 1 1 


1 1 1 1 1 1 






1 




1 1 1 1 1 1 1 1 1 


1 

I i 

1 

II 1 1 111 


1 1 1 1 1 1 


p=5.325 : 

1 1 1 1 1 1 



0.1 - 




200 400 



600 800 
trajectory 



1000 1200 1400 



Figure 28: (IW3I) evol: 16^ x 4, mo = 1.9, = 24, m/ = 0.02 
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Figure 33: {qq) vs. (3: 16^ x 4, mo = 1.9, = 24, m/ = 0.02 
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Figure 34: (11^31) vs. (3: 16^ x 4, mo = 1.9, = 24, m/ = 0.02 
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Figure 35: {qq) vs. T (MeV): 16^ x 4,mo = 1.9, = 24,m/ = 0.02 
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Figure 36: (11^31) vs. T (MeV): 16^ x 4, mo = 1.9, = 24, m/ = 0.02 
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Figure 37: (gg) vs. (3: 16^ x 4, ci=-0.331, mo=1.9, L,=24, m/=0.02 
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Figure 38: (IPF3I) vs. (3: 16^ x 4, Ci=-0.331, mo=1.9, L,=24, m/=0.02 
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Figure 39: {qq) vs. T (MeV): 16^ x 4, Ci=-0.331, mo=1.9, L,=24, m/=0.02 
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Figure 40: {{W^D vs. T (MeV): 16^ x 4, ci=-0.331, mo=1.9, L,=24, m/=0.02 
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Figure 41: Fit of AJ5 = coe-"'^^; 8^ x 4, /3 = 5.2, mo = 1.9, m/ = 0.02 
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Figure 42: Lg dependence of mres and bo'. 8^ x 4, /3 = 5.2, mo — 1.9, m/ = 0.02 



